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Optimal T-S models for identification of nonlinear systems from 

input-output data 

H.Ouakka, I.Boumhidi 

L.E.S.S.I, Department of physics, Faculty of Sciences ,Dhar El mahraz 
BP 1796 FES Atlas, Morocco 
h.ouakka@iam.ma, iboumhidi@hotmail.com 



Abstract 

Determining optimal structure identification of a non- 
linear black box system is one of the most signifi- 
cant steps in Fuzzy modelling based on Takagi-Sugeno 
models. The one parameter that needs to be deter- 
mined before performing fuzzy clustering and iden- 
tification algorithms is the optimal number of clus- 
ters. In this paper, we present a new approach for 
automatically predicting the optimal choice of this pa- 
rameter simply by detecting the general trend of the 
data structure. First, an approximation model of the 
system is built by fitting data to a polynomial func- 
tion. Second, a preliminary decomposition of the data 
is realized based on detection of the function turn- 
ing points. Then, a merging method is adopted to 
reduce the identified number of clusters. The advan- 
tage of the generated solution is that it remains in the 
horizon of the data; hence there is no need to apply 
heuristic rules or conventional validation tools. The 
performance of the proposed method is evaluated for 
both quality of clustering and fuzzy modeling with 1 st 
order TS systems using GK fuzzy algorithm . 

Keywords: Fuzzy clustering, optimal clusters 

number, nonlinear system, polynomial regression, 
Takagi-Sugeno models, GK algorithm. 



1 Introduction 

Fuzzy clustering combining with regression analysis 
are the most tools used to capture sub linear sys- 
tems in a structure of data, particularly for mod- 
elling an unknown black box nonlinear system with 
Takagi-Sugeno models [3] . There are two main is- 
sues in the process of constructing a TS fuzzy model, 
the first is how to determine the number of rules and 
the variables involved in the rule antecedents, and the 
second is how to estimate the parameters of the TS 
fuzzy model. To construct TS fuzzy model from data, 
many fuzzy clustering algorithms combined with least 



square method have proved to be suitable techniques 
to partition data space into adequate subspaces and 
detect linear local models [5] . 

The Fuzzy C-Means(FCM) and the Gustafson- 
Kessel(GK) fuzzy clustering algorithms are the most 
widely used to partition a data set into c homogeneous 
fuzzy clusters. The FCM detects clusters of a roughly 
same size. The GK algorithm is an extension of the 
FCM, which can detect clusters of different orientation 
and shape in a data set . Both the FCM and the GK 
algorithm require the number of clusters as an input, 
and the analysis result can vary greatly depending on 
the value chosen for this variable. However, in many 
cases the exact number of clusters in a data set is not 
known. Both the FCM and the GK algorithm may 
lead to undesired results if a wrong cluster number is 
given. The number of clusters needs to be known a pri- 
ori and the correct number of rules must be estimated 
in advance [1, 10] . 

To overcome this problem, many approaches have 
been proposed: in the first framework, modelling 

procedure constitutes an iterative process, where the 
number of rules gradually increases until the model’s 
performance meets a predefined accuracy[7, 12] , the 
second modelling framework is based on using optimal 
fuzzy clustering, starting with a sufficiently large num- 
ber of clusters and successively reducing this number 
by merging clusters that are similar (compatible) with 
respect to some predefined criteria [11, 13]. 

In this work we present a new technique to deter- 
mine the optimal number of clusters, T-S models, that 
should be used for modelling a given data set of a non 
linear black box system. The basic idea is to fit a 
polynomial function which follows the general trend 
in data. The global structure of data space is then 
recognized by the characteristics points: maxima and 
minima. Each group of data ranging between max- 
imum and minima, considered a cluster, can be ap- 
proximated by a linear model. So the total number 
of cluster, i.e. linear models, can be deduced from 
the number of turning points. To find the optimal 
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number, a cluster merging technique is proposed. It 
consists on fusion of consecutive clusters for witch 
the associated linear model are almost Collinear .The 
performances of results are examined by using the so 
called validity clustering criteria and the quality of the 
approximation is evaluated by comparing fuzzy mod- 
els obtained by the MATLAB Fuzzy Toolbox and the 
fuzzy Model identification based on Gustafson-Kessel 
clustering algorithm. 

The basic advantage of the proposed method when 
compared to other works, is that the optimal num- 
ber of TS models obtained may be used as an imput 
parameter by any TS identification algorithm to gen- 
erate an optimal fuzzy model without need to apply 
heuristic rules or conventional tools to validate the 
results. In more, the peformance of the most indices 
found in the literature depends on the algorithm used 
and the strcutre of the studied system [6, 15]. 

This paper is organized as follows. Structure of TS 
fuzzy model is presented in section 2. Section 3 de- 
scribes the method for construction of the approxi- 
mation polynomial function , the partitioning method 
for extraction of the global number of clusters and the 
merging technique . Section 4 presents simulation re- 
sults . Section 5 gives some concluding remarks and 
discuss current and future extensions of this work. 



2 Takagi-Sugeno fuzzy models 

Consider the identitcation of an unknown nonlinear 
Black Box system of the form : 

Uk = f(%k) 

based on specified or measured input data X and mea- 
sured output data Y of the system, where k = 1, ..., TV, 
denotes the index of the k-th input-output data-pair. 
In general it may not be easy to find a global nonlin- 
ear model that is universally applicable to describe the 
unknown system /(.).In that case it would certainly 
be worthwhile to build local linear models for spe- 
cific region of the data space and combine these into 
a global model. The realtion between input and the 
outupt of the process is then descrribed by linguistic 
if-then rules with fuzzy proposition in the antecedent 
and crisp functions in the consequents. Hence, it can 
be seen as a combination of linguistic and mathe- 
matical regression modeling in the sense that the an- 
tecedents describe fuzzy regions in the input space in 
which consequent functions are valid. The TS rules 
have the following form:: 

Ri : if x is Ai then yi = aix+bi i = l,2,..,c 

where c is the number of if-then rules, Ai are fuzzy 
sets, and (a^bi) are the model consequent parame- 
ters that have to be identified in a given data set. For 
a given input crisp vector x = • • • , £at) t the 

inferred global output of the Takagi-Sugeno model is 




Figure 1: Results of the GK algorithm for fuzzy mod- 
elling of a black box nonlinear system : Fuzzy model, TS 
models and Clusters. 

computed by taking the weighted average of the indi- 
vidual rules contributions 

Ei=i 

y Ei=i ^ 0*0 

where pi(x) is the degree of fulfillment of the i th fuzzy 
rule. 

The Gustafson. Kessel (GK) clustering algorithm has 
often been applied to identify automatically TS mod- 
els. In general , fuzzy clustering algorithms generate 
a fuzzy partition given as a fuzzy partition matrix 
U = [pij], where [lij = fip (xj) is the membership 
value of the data Xj belonging to the fuzzy cluster 
The objective of GK algorithm is to obtain a fuzzy 
c-partition F = {Fi, . . . , F c } for the given number of 
clusters c and the given N input-output data (X,Y) 
by minimizing the evaluation function J m , 

c N 

Jm(U, V,A:X) = Yi 

1=1 j= 1 

D ijAi = ( x j - Vi) T Ai(xj - Vi) 

where V = (rq, . . . , v c ) is a vector of the centroids of 
the fuzzy clusters (Fi, . . . , F c ),m controls the fuzziness 
of membership of each datum , D is an adaptive norm ( 
Mahalanobis distance) used by the GK algorithm for 
each cluster to detect different geometrical shapes in 
data sets. Each ith cluster has its own norm-inducing 
matrix Ai. where Ai is obtained from fuzzy covariance 
matrix (Fi) of the ith cluster defined by: 

p Ef=l ^(Xj-ViUXj-ViF 

' y N u m 

2^j=i M ij 

Ai = [ Pi det(Fi)}- 1/p F- 1 

with pi > 0 Vi, p is the data space dimension . 
Without any prior knowledge, the cluster volume pi 
is simply fixed at 1 for each cluster, which allows 
the algorithm to find clusters of approximately equal 
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volumes [2]. 

It should be noted that for the process of identifica- 
tion of T-S models based on GK fuzzy clustering algo- 
rithm, the number of rules coincides with the number 
of clusters, each cluster is approximated with a local 
linear model(Figurel). 



that can better predict the output value y on the ba- 
sis of the regression vector x is given by : 

V = f(x) = polyfit(x, y, n). (1) 

n is the polynomial order or degree. 

The functioni / is defined by: 



3 Method formulation 

In this section , we present our method in tree steps: 
identification of the polynomial model, extracting of 
the global number of clusters from the approximation 
function and the optimization of T-S models by the 
merging thechnique ( Figure 2) . 



Step 1 



Step 2 



Step 3 




f(x) = a 0 + Y, w k x k (2) 

k= 1 



Depending on the degree n, these models can be 
used to approximate a large class of nonlinear 
functions. However, the number of parameters w in- 
creases with the size of X and the degree n.That 
clearly reflects that polynomial models suffer from 
the curse of dimensionality, problem known as Runge 
phenomenon[4].In general, the estimation of n can be 
performed by stepwise regression. 

To overcome this probleme, an heuristic formula (3), 
is proposed to compute the adequate polynomial or- 
der to build an intermediate mathematical function 
that describes the overall structure of the data and 
preserve features of the distribution such as relative 
maxima, minima and width. 

n = y/N • (1 + r 2 ) (3) 

with r = correlation coefficient (X, Y ) . 

To illustrate this approach, (seefigurel), let us con- 
sider the following process (4) representing a noisy 
SISO black box non linear system : 

y = sinc(x) + e (4) 



Figure 2: Process for identification of the optimal 
number of TS models. 



3.1 Data model identification 



e a vector of disturbance. 

In this exemple ,the value of the polynomial degree 
computed by (3) is n = 12, the model is more powerful 
compared to the others (n = 5 , n = 7), the overall 
aspect of the system is guaranteed and the error of 
modeling is minimized. 



The general problem studied in this section can be 
stated as follows. On the basis of N input-output sam- 
ples 

(xi,yi),(x 2 ,V2), ■ ■ ■ ,(xn,Vn) 



where Xi E X is the ith input vector and yi E Y is 
the corresponding output, find the model / that best 
approximates the true underlying function mapping x 
to y . 

In that way, nonlinear regression by polynomial mod- 
els (using the matlab function polyfit) is used to build 
a nonlinear function approximation that preserve the 
general structure of the data . The regression prob- 
lem stated as a function approximation problem is the 
following. 

From N samples (xi,yi) , the polynomial function 




Figure 3: Nonlinear regression by polynomial models 
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3.2 Extracting clusters from function 
approximation 

Analyzing the behavior of the polynomial function 
graph (figure 3), we can see that the data structure 
can be partitionned in seven substructures or clusters 
.This number is equal to the number of turning points 
, maximum and minimum, increased by unit. 
Consequently, we can propose the following result : 
finding the correct number of cluster involves the 
determination of the number of turning points of the 
approximation function described by the polynomial 
model. 

According to the two fundamentals theorems of 
algebra: 

Theorem 1 Any polynomial P(x) = «o + X^=i w k% k 
of degree n , with real coefficients has precisely n zeros 
in C. 

Theorem 2 A polynomial function of degree n,can 
have at most (n — 1) turning points. 

The total number of turning points are the roots of 
the polynomial function P(x) , determined by using 
the command matlab polyder. In our case, the derived 
polynomial P (x) of degree (n — 1) has (n — 1) real or 
complex roots. 

Assuming that a is the number of real roots and (3 the 
number of complex roots , referring to the first and 
second theorem, the polynomial degree and roots are 
bound by the following equation (5): 

a + j3 = n — 1 (5) 

Since the system is defined in the data space ( X * Y) 

, we are interested thereafter only in the real roots 
which are included in the data input interval X . Sup- 
pose that a* is the number of total real roots A i witch 
define the turning points abscises, such as : 

A = [Ax, ... , A a *] G X 

The data input space can be partitioned in (a*) in- 
tervals: 

X = [x \ , . . . , Ai] U [Ai , . . . , A 2 ] U, . . . , U [A a * , . . . , xn} 

V y v -/ V -/ 

Xi x 2 x c * 

As a result, according to the remark made at the be- 
ginning of the paragraph, the number c* , defined in 
(6), is the global number of cluster that should be 
used to partition the data space (7). 

c* = a* + 1 (6) 

cluster 

X*Y=\jJxfn- (7) 

Xi and Yi are respectively the ith sub input data in- 
terval and its corresponding sub ouput interval. 

By applying the least square method for each group 
of data in a cluster, we can identify for each cluster 
a local linear model: y = a{X + i = 1, . . . , c* that 
approximates each sub system as shown in figure 4. 




Figure 4: Partition of input data set in 7 subsets based 
on function approximation turning points. 

3.3 Determining the number optimal 
of cluster 

In figure 4, the two first regressed linear mod- 
els present a strongly correlation. The corresponding 
groups of data may be merged in one cluster. In that 
way ,for optimizing the number of clusters, we propose 
a clustering merging criterion based on the fusion of 
two consecutive line segment i and i + 1 , given by 
(8) into a single cluster if their associated vectors are 
nearly correlated or linear. We assume that each local 




Figure 5: Rreduction of the number of local models 
(clusters) to 6 by the merging technique. 

regression model is interpreted as a vector Vi defined 
for x G Xi,y G Y im 

y = aiX + bi,i = 1 , . . . ,c*. (8) 

It is known, trigonometry and geometrical approach, 
see figure 4 , that the co-linearity or correlation of two 
vectors is measurable by the angle : 

0 = angle(V i: V i+1 ), i = 1, • • • , (c* - 1) 

Vi is the vector corresponding to the ith local linear 
model. 

In addition there is a relationship between this angle 
and the correlation of the vectors assumed repre- 
senting the regressed data in each cluster cos# = r . 
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The angle for witch the merging criteria is validate, 
is taken equal to 0 r = 30 , this is equivalent for a 
correlation coefficient of data equal to r = 0.82 , 
considered strongly correlation. The cluster merging 
algorithm is presented bellow: 



e is a distribution added to simulate a balck box non- 
linear system. 

The Gustafson-Kessel (GK) algorithm is applied to 
cluster these data set at each cluster number c from 
c = 4 to c = 11 . 



1. Initial parameters 

• k = 0, Fusion number counter. 

• 0 r = 30, V alidation merge criterion . 

for i = 1 to (c* — 1) 

2. Construction of vectors and angle 0. 

• Vi <- y = diX + bi,x G Xi. 

• Vi+ 1 y = \X + bi+ 1, X G Xi+ 1. 

• 0 = angle(Vi,Vi+ 1). 

3. Merge criterion control. 




• if 0 <6 r . 

• Merging intervalsX i = Xi U X i+ i . 

• Resulting local linear model y = a^Fb^ 

• k = k + 1. 

The optimal number of clusters corrsponding to 
the optimal partition is then given by: 

C = c* -k. (9) 

The number of cluster determined ,in subsection 
3.2, for the system (4), is equal to seven clusters; 
this number is reduced to six by applying the merging 
method. Figure 5 shows the new partitionnig of the 
X data input ,the sub data inputs X\ , X 2 are merged 
and a new sub local linear model is built to approx- 
imate the data in the resulting cluster corresponding 
to sub input interval [xi, , A2]: 

X = [xi , . . . , A 2 ] U [A 2 , . . . , A 3 ] U, . . . , U [A6, • • • , xn] 

X[ X' X' 

X[ is the resulting sub input interval from the fusion 
of X\ and X 2 . 

Comparing the fuzzy model ( fgure 1) and the poly- 
nomial function ( figure 3) , the two models represent 
identically the general trend of data, this shows well, 
the effectiveness of the polynomial function to be used 
as a basis for partitioning of the data of an unknown 
nonlinear system. 

4 Simulations results 

In order to verify the results of the proposed 
method, we have prepared a set of noisy data gener- 
ated by the function (10) taken from [9, 8] approxi- 
mated from 201 points placed in equal intervals in [-1, 
!]• 

y = 0.6 sin(7nr) + 0.3 sin(37nr) + 0.1 sin(57nr) + e. (10) 



Figure 6: Result of partition of the system (10) data 
set : polynomial model ( n=ll ), Optimal number of 
clusters (c=7). 



Tow validations approaches will be used to eval- 
uate the results. The quality of clustering is tested 
using the fuzzy clustering toolbox to compute the so 
called clustering validity criteria : Partition coefficient 
(PC) and Classification entropy (CE) .Both the two 
indices measures the amount of overlapping between 
clusters. The optimal partition can be determined by 
the point of the extrema of the validation indexes 
in dependence of the number of clusters. The optimal 
number of clusters is assumed equal to the local max- 
imum of the (PC) indices and the local minimum for 
the (CE)[14]. 

The performance of the fuzzy model is measured by 
the VAF, which computes the percentile Variance ac- 
counted for between two signals as follows: 



VAF = 100 • 



varjyi - m) 
var(yi) 



7/1 is the output of the process and y 2 is the output 
of the fuzzy model with the same number of clusters 
obtained by the MATLAB fuzzy toolbox ( ANFIS ) 
and the fuzzy model identification ( FMID) toolbox 
based on GK clustering[l]. The VAF of two equal 
signals is 100. If the signals differ, VAF is lower. 

Figure 6 illustrate the result of our proposed method 
for partitionnig the data based on the identified 
polynomial function .The plynomial degree computed 
by (3) is n = 11 , number of turnings points detected 
is 6 and the optimal number of clusters is given by 
(9) c = 7. 

The performance of each clustering validity measure 
is given in figure 7. As seen from the figure, the (PC) 
and (CE) validity measures find that the optimal 
cluster number c is at c = 7 . 
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4 5 6 7 8 9 10 11 12 

Number of clusters 



Figure 7: Performance of (PC) and (CE) clustering 
validation index for different number of clusters. 




Number of clusters 



Figure 8: VAF results for different number of clusters. 

In figure 8, The VAF predict c = 7 as the correct 
optimal number of clusters since there is no wide 
variation between the value of the VAF for the 
number of cluster ranging between 7 and 11. This 
shows applicability of our method to compute the 
optimal number of clusters for approximation a 
nonlinear black box system with a fuzzy model as 
shown in figure 9. 



nonlinear system, the system (10) is taken without 
perturbation, our method has computed (c = 9) opti- 
mal clusters ,the algorithm using a complementary in- 
terpolation model applied in [9] has validate the same 
number of interpretable TS models that capture the 
essence of the process under consideration. 

The simulation results are illustrated in table 1. 



cluster 

number 


VAF 

FMID 


VAF 

ANFIS 


PC 


CE 


6 


98.208 


99.200 


0.840 


0.309 


7 


99.613 


99.912 


0.868 


0.264 


8 


99.705 


99.941 


0.881 


0.234 


9 


99.821 


99.967 


0.894 


0.207 


10 


99.760 


99.961 


0.892 


0.212 


11 


99.974 


99.972 


0.894 


0.207 



Table 1: Validation indicies for the system (10) with- 
out pertubation . 



5 Conclusion 

The major contribution of this paper is an improved 
method to compute automatically the optimal num- 
ber of clusters, needed as input parameter for cluster- 
ing algorithms, for correctly approximating a nonlin- 
ear black box system with TS models. The simulation 
show that the proposed approach simplifies and makes 
more effective to extract optimal and legitimate TS lo- 
cal models from data . In addition, the data model 
identification used in this paper has proved to be suit- 
able for transforming an unkown balck box system 
into a grey box one . The approach can be extended 
for more complexes systems, since it has proved im- 
portant result for SISO nonlinear systems. 
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Abstract 

This article presents sufficient conditions for the 
uniform stability of the zero solution when p = 0, 
the uniform boundedness and uniform ultimate 
boundedness of all solutions in case p fz 0 of the 
third-order nonlinear delay differential equation (1). 
By constructing Lyapunov functionals we achieved 
new contribution which includes and improves some 
related results existing in the relevant literature. 

Keywords: Lyapunov functional, Stability, 

Boundedness, Third- order delay differential equa- 
tions. 



1 Introduction 

Delay differential equations was encountered in the 
late eighteenth century by the Bernoulli’s, Laplace 
and Condorcet. However, little was accomplished dur- 
ing the nineteenth century and early part of the twen- 
tieth century. During the last forty years and espe- 
cially the last twenty, they have applications in do- 
mains as diverse as engineering, biology and medicine, 
where information transmission and/or response in 
control systems is not instantaneous. 

Stability is a very important problem in the theory 
and applications of differential equations. The stabil- 
ity analysis of delay differential equations has received 
considerable attention over the decades. The stabil- 
ity of delay differential equations was considered in 
[7, 16]. Also boundedness of solutions was discussed 
in [16]. Later many books and papers dealt with the 
delay differential equations and obtained many good 
results, for example, [1, 3, 4, 5, 6, 8, 10], etc. 

Lyapunov function is an interesting and fruitful 



technique to determine the stability behavior of so- 
lutions of linear and nonlinear differential equations. 
This technique has gained increasing significance and 
has given impetus for modern development of stability 
theory of differential equations. 

Besides, it is worth mentioning that, according our 
observation there are only a few papers on the same 
topic for certain third-order ordinary nonlinear dif- 
ferential equations with delay. Perhaps, the possible 
difficulties raising this case are due to construction of 
Lyapunov functional for the delay differential equa- 
tions. 

In fact, there is no general method to construct the 
Lyapunov functions and clearly, it is more difficult to 
construct Laypunov functionals for higher-order dif- 
ferential equations with delay. 

To the best of our knowledge, for a kind of third- 
order nonlinear delay differential equations, the sta- 
bility and boundedness results have been investigated 
only by a few researchers, for example, In 2003, Sadek 
in] obtained sufficient conditions to ensure the stabil- 
ity and the boundedness of solutions of the equation 

a; (3) +ax + g(x(t - r{t ))) + f(x(t - r(t ))) = p(t), 

where a is positive constant; g,f,p are continuous 
functions; g(0) = /( 0) = 0. 

In 2004, Sadek [12] investigated the asymptotic 
stability of the zero solution of the fourth-order non- 
linear delay differential equations 

x^ + ol{x + a2X + asx + f(x(t — r)) = 0, 

X^ + OL\X + OL^X + cj)(x(t — r)) + f(x) = 0, 

where r, aq, ot2, ol 3 are positive constants; (f>{x ) and 
f{x) are continuous functions and 0(0) = /( 0) = 0. 

In 2005, Sadek [13] gave the sufficient conditions 
for the asymptotic stability of zero solution of third- 
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order delay differential equation 



2 Stability 



x ^ + a(t)x + b(t)x + c(t)f(x(t — r)) = 0, 

where a(£), b(t) and c(t) are positive and continuously 
differentiable functions on [0, oo); r is positive con- 
stant; f(x) is continuous function and /( 0) = 0. 

In 2006, Tung [14] obtained the asymptotic stabil- 
ity of solutions of third-order delay differential equa- 
tion of the form 

x x)x + g(x(t - r(t))) + f(x(t - r(t))) = 0, 

and the boundedness of the equation 

x + ax + g(x(t — r(t ))) + f(x(t — r(t ))) = p(t, x, x, x) 

where 0(x, x), g(x), f(x) and p(t,x,x,x) are continu- 
ous functions; g(0) = /( 0) = 0. 

In 2007, Tung [15] investigated stability and 
boundedness of solutions of some nonlinear differen- 
tial equations of third-order with delay of the type 

x (t)+aix (t) + f 2 (x (t - r(t))) + a 3 x(t) 

= p(t, x(t),x(t),x(t - r(t)),x(t - r(t)), x (£)), 



By constructing Lyapunov functionals, we obtained 
very clear results which can be checked and applied 
easily. The following is the first main result. 

Theorem 1. Let h(0) = g(x, 0) = 0 for all 
x. Suppose that there exist positive constants 
a, 6, c, L, M, N and ao with ab — c > 0 such that 

(*) 0 < f(x, y)-a<a 0 , y 9/ ^’ y) < 0 for all x, y, 
(ii) yy >b> 0 ,y^ 0 -, 

(Hi) sup {h'(x)} = c > 0, h(x)sgnx >0 for x ^ 0 ; 

( iv ) r(t) < 7, r'(t) < /?, 0 < /3 < 1; 

(v) \h’(x)\ < L, I |< M, I |< AT. 

Then the zero solution of (1) with p = 0 is uniformly 
stable, provided that 

ry < min\ (afo— c— 2M)(1— /3) 

7 <. nun /x ( L+M+A r)(i_ / 3) + ( L+M )(i +M ), 

(ab-c)(l-P) 

6{(L+M)(l-/3)+AT(/x-/3+2)} ’ 



where r is abounded delay, 0 < r(t) < 7, r(t) < /?, 
0 < /3 < 1, /3 and 7 are some positive constants; 
a\ and as are some positive constants; f 2 and p are 
continuous functions and /2(C)) = 0. 

Here we consider the third-order nonlinear delay 
differential equation of the form: 

x ( 3 ) _j_ ±)x + g(x(t — r(t)),x(t — r(t))) 

+h(x(t — r(t ))) = p(t, x, x, x, x(t — r(t)),x(t — r(£))), 

(i) 

where 0 < r(t) < 7, 7 is positive constant which will 
be determined later; f,g,h and p are continuous func- 
tions; h( 0) = g(x, 0) = 0. 

Our object is to obtain sufficient conditions for the 
stability and for the boundedness of solutions of (1) 
in the cases p = 0 and p 7^ 0 respectively. 

The results of this paper generalize those obtained 
in Sadek [11], where f(x,x) is a constant, 
g(x(t — r(t)),x(t — r(t ))) was considered as a function 
of x(t—r(t)) only andp(£, x, x , x, x(t—r(t)),x(t—r(t))) 
was considered as a function of t only. 

The remainder of this paper is organized as fol- 
lows: 

Section (2) gives the sufficient conditions for the sta- 
bility of the zero solution of equation (1). 

Section (3) shows the boundedness of all solutions of 
equation (1). 

Section (4) presents example which illustrate the de- 
sign of nonlinear controllers in feedback systems via 
Lyapunov’s criterion. Proceeding from the knowledge 
of Lyapunov functions and their stability conditions 
the form of the nonlinearities is obtained and their ef- 
fect on the overall nonlinear feedback system is stud- 
ied. 



where fi = (ab + c)/ 2 b. 

Now we will give the stability criteria for the general 
autonomous delay differential system. We consider: 

x = f(xt), x t = x(t + 0), —r< 6 < 0, t > 0, (2) 

where / : Ch — > is a continuous mapping, 

/( 0) = 0, C H := {(/> e (C[-r, 0],5R") : ||0|| < H} and 
for Hi < H , there exists L(H\) > 0, with 1/(0) I < 
L(Hi) when ||0|| < H x . 

Lemma l.[l, 17] Let V(4>) : Ch — ► 3? be a contin- 
uous functional satisfying a local Lipschitz condition, 
V/O) =0 and functions H/(r), (i = 1,2) are wedges 
such that 

(*) W^mi) <V{ 4 >) <W 2 m\) and 

(ii) V( 2 )( 4 >) < 0, for e C H - 

Then the zero solution of (2) is uniformly stable. 
Proof of Theorem 1. 

We write equation (1) with p = 0 as the following 
equivalent system: 



x = y, 
y 

z ■ 



~f(x,y)z - g(x,y) - h(x) 

+s; r m 

+ ft-rlt) h’(x(s))y(s)ds. 



Define the Lyapunov functional as 

Vi (x t ,y t ,z t ) = n fg h(£)d£ + h(x)y 

+M fo f( x ’ V)vdy + Jo g(x, rj)drj + yyz 



+h 2 + A/7 (t) / t 

+ 6 j°_ r{t) j t t+s z 2 ( 0 )deds, 



r{t)Jt+ s y 2 (0)deds 



(3) 



(4) 
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where fi, A and S are positive constants which will be 
determined later. Thus from (4) and (3) we find 

$tVi(x u yuZt) = h'(x)y 2 + yz 2 - yyg(x,y) 

-f(x, y)z 2 + y fo y ^^- ydri + ff y ^-f^ dy 
+{yy + 2 ) | / t *_ r(t) h'(x(s))y(s)ds 

+/.-.(,) a,( y » »(■)<<»} 

+\y 2 r(t) - A(1 - r'(t)) //_ r(t) y 2 (9)d9 
+Sz 2 r(t) - <5(1 - r'{t)) / t *_ r(t) z 2 (9)d9. 

(5) 

From (i) and (v) we find 

j- t Vi ( x t,yt,z t ) < h'{x)y 2 + yz 2 - yyg(x,y) + My 2 
-f(x,y)z 2 + (yy + z) j(M + L) //_ r(t) y(s)ds 

+N z{s)ds I 

+Ay 2 r(t) - A(1 - r'(t)) / t *_ r(t) y 2 (9)d9 
+5z 2 r(t) - (5(1 - r'{t)) / t *_ r(t) z 2 (9)d9. 

(6) 

From (i), (ii), (Hi), (iv) and by using 2 uv < u 2 + v 2 
we obtain 



& U < ^ L+ f +JV) 7 -A 7 |y 2 

— ja — y — L+ ^ +jV 7 - £7 jz 2 

+ j L±M + M^+M) _ (1 _ /3)A j y 2 (s)ds 

+ |t + £ T' _ ( 1_ j ft—r(t) Z 2 (s)ds. 

(7) 

We can take y = > o, A = > 0, 

5 = > 0 so that 



^L+M+AQg-ffl + fL+MHp+l) ) 2 

2(i-/3) 

_/ a^c _ (L+M)(l-/3)+iV( M -/3+2) .,t „ 2 
| 2b 2(1— (3) I p • 



Therefore if 



7 < min 



(ab—c—2M)(l—/3) 
M(L+M+AT)(l-/3) + (L+M)(l+/i) ’ 



(afr-c)(l-/3) 

6{(L+M)(1-/3)+AT( m -/3+2)} > 



From (in) by considering h(x)sgnx > 0 for x / 0, we 
obtain 



Vi(x t , y t , zt ) > y ff h(£)d£ + h(x)y + \yay 



-\by 2 + yyz + ±z 2 + A f° f 



+<5 f°_ 



f t , s z 2 (0)d6ds 



r(t) Jt+s 



y 2 (6)d0ds 



r(t) Ji+s 



= ib( b y + h ( x )) 2 + \{yy + z) 2 + \y(a - y)y 2 

4 fo h (0\ Jo (y b ~ h\i))ydrXdi 



2b 

+ 



2 by 2 



+xj: r{t) j; +3 y 2 (9)d9ds 

+<5 f-r(t)Jt+s z2 (°) d0ds ’ for y^°- 

But a — fi = (ab — c)/2b > 0, 

gb — h'(x) > (ab — c )/ 2 > 0 and since the integarls 
A fl r(t) Jt +s y 2 (9)d0ds and (5 /° p(t) f* +s z 2 (9)d9ds 
are non negative, then we obtain 



Vi(xt,yt,z t ) > £( by + h(x )) 2 + \(yy + zf 
+\y(a - y)y 2 , for y ^ 0. 



(9) 



From (i), (iv), (v) and by using the mean- value theo- 
rem we find 



Vi (x t , y t , z t ) < y fo L£d£ + L\xy \ + y ff (a 0 + a)ydg 

+ fo Nr i d v + 7 2 + mM 

+A ft_ r(t) (9 -t + r(t))y 2 (9)d9 
+(5 fl_ r{t) (9 -t + r(t))z 2 (9)d9 
< (M+ 2 1)L bH 2 + L+JV+ ^( a + a o + 1 ) + W 2 || y ||2 
+B± 1±S£\\ Z \\ 2 ' 

( 10 ) 

Then from (8), (9) and (10) V\(x t ,yt,zt) satisfies the 
conditions of Lemma 1. 

Thus the proof of Theorem 1 is now complete. 



3 Boundedness 

The following is the second main result. 

Theorem 2. Suppose, further to the conditions of 
Theorem 1, that there exist constants Co > 0, m > 0 
and mo > 0 such that > Co for 
p < m + mo(\x\ + \y\ + \z\). Then every solution of 
equation (1) is uniformly bounded and uniformly ul- 
timately bounded, provided that, 

7 <mmj L+ M> +JV , 

(ab—c—2M)(l—/3) 

(L+M)(l+/x+a6-c+a+a 2 ) + (L+M+AT)( / i+a 2 )(l-/3)’ 

(ab—c)(l—f3) 1 

6AT(l+/i+a6-c+a+a 2 )+6(L+M+AT)(l+a)(l-/3) f ’ 

where fi = (ab + c)/ 2b. 

Now we will consider a system of delay differential 
equations 



we have 



x = F(t, x t ), x t = x(t + 0), —r<6< 0, (11) 



j- f Vi{x t ,y t ,z t ) < -a(y 2 + z 2 ), for some 



a > 0. where F : 5ft x C — > 5R n is a continuous mapping, and 
(8) takes bounded sets into bounded sets. 
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The following Lemma is a well-known result obtained 
by Burton [2]. 

Lemma 1. Let V(t,(j)) : 5ft x C — > 5ft be continuous 
and locally Lipschitz in </>. If 

(*) W(\x(t)\) < V(t,x t ) < Wrdxm + 

W 2 ^Jl r W 3 (\x(s)\)d s y and 



(ii) V(n < —Ws(\x(t)\) + M, for some M > 0 
where W(r), Wi(r),(i = 1,2,3) are wedges. 
Then the solutions of (11) are uniformly 
bounded and uniformly ultimately bounded for 
bound B. 

Proof of Theorem 2. 

Now we will consider the boundedness of the solutions 
of equation (1) by using this lemma. 

We assume that p(t, x, x, x, x(t — r(t)),x(t — r(t))) is 
bounded with bound m + mo(\x\ + \y\ + \z\) and the 
conditions of Theorem 1 hold. Equation (1) has an 
equivalent system 



x = y, 

V = z, 

z = -f 0, y)z - g(x, y ) - h(x) 

+ Sl m djM ^ y(s)ds 

+ ltr(t) 9 -^§jf^ z(s)ds 

+ ft- r (t) h'(x(s))y(s)ds 

+p(t, x, y, z, x(t - r(t)), y(t - r(t))). 



( 12 ) 



Consider the Lyapunov functional 



V = Vi (x t ,y t ,z t ) + V 2 (x t ,y t ,z t ), (13) 



where Vi(x t ,y t ,z t ) is defined as (4) and V 2 (x t ,yt,z t ) 
is defined as 



V 2 (x t , y t , zt ) = a 2 f* h(£)d£ + a 2 / Q y f(x, r^ydy 
+| (ab - c)x 2 + ah(x)y + § y 2 
+ (ab — c)x(z + ay) + | z 2 + a 2 zy. 

(14) 

From (4), (7) and (12) we find 



< — | ^ - M - At(L +f+ 7V) 7 - A 7 jy 2 

_|a|-c _ L+M+N^ _ s ^j z 2 

+ j L±M + _ (1 _ £)* J y 2 (s)ds 

+ |T + i T'-( 1- ft-r(t) Z 2 (s)ds 

+(m|2/| + \z\){m + m 0 (\x\ + |y| + N)}. 



(15) 



Also from (14), (12) and by using the conditions in 

Theorem 1, we obtain 

jr t V 2 (x t ,y u z t ) < -(ab - c)xh(x) 

+(ab - c)\x\{m + m 0 (|x| + \y\ + |z|)} 

+a\z + ay\{m + m 0 (\x\ + \y\ + |z|)} 

+(ab - c)xN f*_ r(t) z(s)ds ( 16 ) 

+(ab - c)x (M + L) f t _ r{t) y(s)ds 
+a(z + ay)N z(s)ds 

+a(z + ay)(M + L) y(s)ds. 

Since r(t) < 7, and by using 2 uv < u 2 + v 2 we have 

jjV 2 (x t ,y t ,z t ) < -(ab - c)xh(x) 

+(ab - c)\x\{m + m 0 (|a;| + \y\ + |z|)} 

+a\z + ay\{m + m 0 (\x\ + \y\ + |z|)} 

+ 2 7 £ (L + M + N)^x 2 + ^(L + M + N)^y 2 

+ f(L + M + N)-yz 2 + | ^(L + M) 

+ f (L + M) + \(L + A4")| St- r (t) y 2 (. s )ds 

+ + f N+ ^n) Jl r(t) z 2 (s)ds. 

(17) 

Therefore from (13), (15) and (17), we obtain 



&V<- 



— M — \~f — 






7 ?y 



J ab—c 

l 2 

-| S^c- L±M±jV 7 - fry - < L+ f f7V) 7 |x 

— (ab - c)xh(x) + <^p(L + M + N)^x 2 
+ | 2*=c(L + M) + ±1(L + M) 






~ (1 ~/5)a| St_ r(t) y 2 (s)ds 

+ | s^N + | 

«(1 -/?)<*} z 2 (s)ds 



+(ab — c)\x\(m + mo) 

+ (a\z + ay \ + y\y\ + \z\)(m + m 0 ). 

Suppose that > Co for x 7^ 0. Then 

ji 



(18) 



f t v < -|(o6 -c)(c- ^±M±AL 7 ) 

|a + ^(1/ + M + AT)(/x + a 2 ) J- 7 



ahp-M 



ab—c 

2b 



- 1 5 + \(L + M + N)(l + a) j 7 J 
(1 — /3) A — l+ 2 m (1 + /u + a6 — c + a + a 2 ) j> 
Arp) y 2 (s)ds- |(l-/3)<5- f(l+/r + a6-c 
+a + a 2 ) | A r(t) z 2 (s)ds + (ab - c)m\x\ 



+ (a 2 m + fim)\y\ + (am + m)\z\ 
+(ab - c)\x\m 0 (\x\ + \y\ + \z\) 

+ (a 2 + y)\y\mo(\x\ + |y| + \z\) 
+ (a + l)| 2 )|mo(|x| + \y\ + |z|). 
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Let A — 2{i— v3) (1 T y H - ab — c a - f- a 2 ) > 0, 

S = 2 (i-p) (1 + /i + a6 — c + a + a 2 ) > 0. When 

7 < min j L+M ° +Ar , 

(ab— c— 2M)(1— /3) 

(L+A4) (l+/x+ab— c+a+a 2 )-|-(L+A4+iV) (/x+a 2 ) (1— /3) 5 

(ab— c)(l— /3) \ 

bAr(l+ A i+ab-c+a+a 2 )+b(L+M+AT)(l+a)(l-/3) J ’ 

where y = (aft + c)/2b. 

Then we can take 

k = mo max (ab — c, a 2 + y, a + l) 

and 

ft (M + 1 2/| + |^|) 2 < 3 k ( x 2 + y 2 + z 2 ) 

We get 

< -a(x 2 + y 2 + z 2 ) + fca(|a;| + |v| + |z|) 

= -f (x 2 + y 2 + £ 2 ) 

-f{(N-fc ) 2 + (l!/l-fc ) 2 + (N^fc) 2 } + ¥^ 

< — f (x 2 + y 2 + z 2 ) + ^/c 2 , for some a, k > 0. 



A direct feedback control system is shown in figure 1. 
( see Ku, Mekel and Su [9]). The system to be con- 
trolled is denoted by its transfer function G\(s). The 
controller is denoted by N in parallel with an ampli- 
fier which gives gain K\. The controller may contain 
linear as well as nonlinear functions of the variable 
x and its derivatives. In an autonomous system, the 
reference input is zero (r = 0). Hence the output C 
is equal to —x, where x denotes error or actuating 
signal. The input to G\(s) is denoted by m, the ma- 
nipulated signal. As shown m is the sum of K\x and 
the output of the controller N, which may be denoted 
by n(t). 

Let 

G l( S ) = sO^+cL+aa)' ( 21 ) 

The amplifier has a gain K\. The differential equation 
of the closed-loop linear system is given by 

x + aq x + a 2 x + K\ x(t — r(t)) = 0, (22) 

controllers with delay. 

Let the nonlinear controller give its output as 

n{t) = fi(x,x) x + gi(x) + h x (x(t “ r(t))). (23) 



Then condition (ii) of Lemma 2 is satisfied by taking 
W 3 (r) = fr 2 and M = ^k 2 . 

From (i) and (in) we have 

V 2 (x t ,y t , z t ) > a 2 f* h(£)d£ + ah(x)y + \o?y 2 
+ |(a& — c)x 2 + %y 2 + (ab — c)x(z + ay) 

+ | z 2 + a 2 zy 

4 /o x Me{/ 0 y (c-/i'(0)^}^ 

+b( c V + ah(x )) 2 + ^(z + ay ) 2 
+ s ^ £ {bx + (z + ay)} 2 . 

(19) 

By using the similar techniques of (10) we obtain 




V(x t ,Vt,Zt) < | (J!±2±^± m + (ab-c)(b+l) 1^2 

J (d— |— do ) (a+l)L+/x-|-iV+a 2 (q-|-q 2 )(qb — c) 



\ 



I / /i+a+q 2 + l (a+l)(ab— c) j ^2 

+^{ //-rw f ON 2 + IMI 2 + IMI 2 ) d «}- 



(20) 



where /3 = max(A, S). 

Therefore from (9), (19) and (20) the functional 
V(xt,yt,zt) satisfies condition (i) of Lemma 2. 

Thus the proof of Theorem 2 is now complete. 



4 Case of Study 

Lyapunov methods have been used by control system 
designers to obtain stabilizing feedback controllers for 
nonlinear system. 

We apply Lyapunov’s direct method to the stabil- 
ity study of nonlinear feedback control systems. 



Then 

m(t) = Ki x(t - r(t)) + fi(x, x) x + gx(x) ( . 

+h 1 (x(t-r(t))). ^ > 




Figure 1: Block diagram of nonlinear system 

The differential equation for the overall nonlinear con- 
trol system with direct feedback is 

£ (3) + (ax + fx(x, x)) x + a 2 x + gx(x) , . 

+hx(x(t - r(t))) + Ki x(t - r(t)) m 0. ^ ’ 

If we let 

f(x,x) = ol\ + fi(x,x), g(x,x) = a 2 x + gi(x), 
h(x(t - r(t))) = K\ x(t - r(t)) + hi (x(t - r(t))), 
we obtain the following nonlinear delay differential 
equation of the third-order 

x ( 3 ) ±) x + g(x, x) + h(x(t — r(t))) = 0. 

(26) 

We note that equation (26) is a special case of equa- 
tion (1) and the functions f,g,h satisfies the condi- 
tions of Theorem 1, provided that <7i(0) = 0 and 
9 -^>b,y^ 0. 
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5 Conclusion 

The most effective method to determine the stability 
behavior and boundedness of solutions of linear and 
non-linear differential equations is the Lyapunov’s di- 
rect method. The major advantage of this method is 
that stability and boundedness in the large can be ob- 
tained without any prior knowledge of solutions. To- 
day, this method is widely recognized as an excellent 
tool not only in the study of differential equations but 
also in the theory of control systems, dynamical sys- 
tems, systems with time lag, power system analysis, 
time varying non-linear feedback systems, and so on. 
This paper gives an introductory study on the sta- 
bility, boundedness and design of nonlinear control 
systems of third-order delay differential equation. It 
shows that Lyapunov’s criterion can be applied with 
advantage from the control engineer’s point of view to 
nonlinear feedback systems. 
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Abstract 

This paper addresses the problem of designing a fuzzy 
logic PI controller for a class of Multi-input Multi-output 
(MIMO) systems having a polynomial input nonlinearity. 
The MIMO system is decomposed into Multi-input 
Single-output (MISO) subsystems. By this 
decomposition, the identification, decoupling and control 
become more performing and easy. We propose to 
approximate ah the interactions between the loops, the 
nonlinearities and disturbances using the polynomial 
approach. The control is conducted at the level of each 
subsystem and the decoupling is easily achieved by 
compensation of the cross-coupling. Then the control of a 
MIMO system is transformed into the control of a set of 
coupled nonlinear MISO subsystems. A good 
performance for tracking is obtained using fuzzy logic PI 
controller. The validity and robustness of the proposed 
approach is tested on a simulation example. 

Keyword: MIMO systems, Adaptive control, Nonlinear 
systems, Decoupling, Fuzzy logic controller 

1. Introduction 

The majority of process industries are nonlinear and 
multivariable. Due the interconnections, dead time and 
nonlinearity, the control of these systems become 
difficult. It is obvious that the difficulty of MIMO 
systems control is how to overcome the coupling effects 
among each degree of freedom. To obtain good 
performance, coupling effect cannot be neglected. Hence 
SISO system control scheme is not easy to implement on 
complicated MIMO systems. Therefore, intelligent 
control strategy is gradually drawing attention. A 
multivariable decoupling control strategy should be 
adapted to obtain good performance for tracking purpose. 
A multivariable adaptive control ensuring decoupling has 
been the topic of several recent papers. On approach is a 
parameterization method to avoid the interactor matrix 
[1,2]. Another approach is to use a precompensator in 
cascade with the plant so that we obtain a diagonal 
interactor. Thus, the existence of such precompensator is 
characterized by a priori knowledge on the relative 
degree of each element of the plant, while the parameter 
values are unknown. However this parameterization is 



quite complicated. Different decoupling approaches can 
be found in the literatures [1, 3, 4, 5, 6, 22]. 

However, ah these approaches are computationally 
intensive and can not be applied to MIMO nonlinear 
systems subject to unmodelled dynamics. 

Fuzzy logic control methodologies have emerged in 
recent years as a promising. Fuzzy logic control has 
found extensive applications for plants that are complex 
and ill-defined. Based on the universal approximation 
theorem [17] several stable adaptive fuzzy control 
schemes have developed to incorporate the expert 
knowledge systematically for a class of nonlinear 
systems with unknown dynamics and disturbances [23, 
24, 25]. The MIMO systems usually process 

characteristics of nonlinear dynamics coupling [2, 5, 11, 
12] and the difficulty is how to overcome the coupling 
effects among each degree of freedom and assuring 
tracking performance. Due to the above problems, we are 
motivated to design a fuzzy logic controller and an 
auxiliary control to compensate the interconnections and 
nonlinearity terms. 

In this paper, the multivariable control design is 
transformed into a monovariable one by considering the 
MIMO system as composed of MISO subsystems. Each 
MISO subsystem will be controlled independently of the 
others. The advantage of this approach is beneficial in the 
sense that it allows the direct practical implementation 
for a large class of multivariable control systems and 
avoids the burden computational of the global approach. 
A tuned fuzzy logic controller FEC is designed according 
to optimize the parameters of the fuzzy system. To 
predict ah the interactions, the nonlinearities and 
disturbances, we propose to use the Chebyshev 
polynomials approach. This prediction is taken into 
account in the local control design and allows to convert 
the nonlinear problem into a linear one. We note that 
discrete orthogonal polynomials have been applied for 
system analysis, parameter identification and optimal 
control with some success [7,8,9]. 

The paper is organized as follows: section 2 describes the 
problem statement. Estimation parameter algorithm is 
given in section 3. The structure of the fuzzy logic PI 
controller is detailed in section 4. In section 5, simulation 
results are given to test the validity of the approach. 
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2. Problem formulation 

We consider a MIMO discrete time system having a 
polynomial input nonlinearity represented by the 
following difference equation [10]: 

A(q- 1 )y(k) = B d (q 1 )Z(k) + m (1) 

Where 

q~ l is the back shift operator : q~ l (x(/:)) = x(k — 1) . 
A(q 1 ) and B d (q 1 ) are polynomials matrix defined 
by: 

A(q~') = diag[A t (q~')\ 

n a (0 (2) 

A M l ) = Yj a i r q r ^ a io= l 

r = 0 

B d (q- l ) = [q- d *B ij (q- 1 )] 

n b (i,j ) (3) 

s„(<r')= 2>,„<r. t> v *o 

r=0 

d u is the time delay between the input u i (k) and the 
output y • ( k ) . 

d- is the time delay between the input U-{k) and the 
output y • (fc) . 

y(£) G , Z(k) G 7?" and g(k) G 7T are 

respectively the output, nonlinear input and the 
disturbances vector. 

The non linearity is assumed to be represented as follows 

[ 11 ]: 

Z ; (0 = fto + fn u i ( k ) + - + f ip u r' ( 4 ) 

The MIMO system (1) can be decomposed into n MISO 
subsystems. The identification, decoupling and control 
will be conducted at the level of each subsystem. Thus, 
the computational load is reduced. 

Equation (1) becomes: 

A i (q- l )y i (k) = q- d “B ii (q- l )Z i (k) + 



A (<r‘ )>’, (k) = q d “ f n B u (q l )u i (k) + 

q~ d “ B n + fn u f (k) +■■■+ (£)] ( 6 ) 

+ f j q^B ij (q-')Z j (k) + ^ i (k) 

7=1 

The objective is to cause the system output 

track a bounded reference sequence y jm (k) using the 

fuzzy logic PI controller. However, it is hardly ensured 
when the system is nonlinear multivariable with 
disturbances. 

The system (6) can be written as: 

A (q-'hi (k) = q- dii b: (q~ l X (k) + H i (.) (7) 

Where 

B*i (q~ l ) = fn B u (q 1 ) (8) 

H • (.) is a non linear function of U i (k ) and Uj(k). 

HA) = q~ dii f i0 B i M~ l ) + q- d ‘ i B ii (q- 1 )Zf ik A(k) + 

k=2 

(9) 

^q^ByiqXZjiV + ^k) 

7=1 

m 

We propose that H i (.) will be predicted in real time 

using Chebyshev orthogonal polynomials. This 
prediction converts the nonlinear problem into a linear 
one. 

In general, the time function f t (t) which is square 
integrable in the time interval 0 can be 

expanded into the orthogonal polynomials. In this paper, 
we use the Chebyshev polynomials to predict /■ (t) in 

real time as: 

00 

/,(»)- Z/.HW (10) 

k = 0 



f d q- d 'By(q- 1 yZ J (k) + Z l (k) (5) 

7=1 

j*i 

Where u i (k) and y • (k) are respectively the i rh scalar 



input and scalar output of the system at time instant k . 
The proposed control law is composed of two terms. The 
first is a PI regulator synthesized using a fuzzy logic 
technique witch allows a good tracking. The second term 
is a compensator, where the structure is obtained from the 

result of the predictor of the function//^.) defined by 



(9). 

Thus, equation (5) becomes: 



which contains an infinite number of terms. In practical 
applications, only a finite number should be used. Thus, 
approximation gives: 

m 

/,(0 = £/,«7-,(0 + cr,m 

k=0 (11) 



where (7 i ( t ) is the truncation error. 



The unknown function //.(.) contains the nonlinear 
terms, cross-coupling and disturbances. Here, we propose 
to estimate H i (.) in real time by projecting it on a set of 
Chebyshev orthogonal polynomials. The estimate of 
7/ ? .(.) will be used in the control design in order to 
synthesise a compensator. 

Based on Chebyshev approximation of //.(.), the system 



described by equation (7) becomes: 
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4 (9 -1 )?,■(£) = q d “B* i (q~ 1 )u i (k) + 

hjTi^ + a^k) 



(12) 



Where CT- ( k ) is the truncation error. 

Remark: concerning the expansion order of the 

Chebyshev polynomials. 

From equation (11), the question which naturally arises is 
how should the expansion order m be chosen? If the 
problem was in a framework of signal analysis, the 
expansion order m should be large enough for good 
signal approximation (theoretically m — » oo ). In the 
present context (from the control point of view), this 
order is chosen to be minimal m = 3 . This is because 
equation (11) is used as a prediction model. It is known 
that augmenting the order m to m + q gives a 
parameter vector (equation (13)) with q extra 
parameters to be estimated, which increases the 
convergence time. 

The advantage of Chebyshev approximation is that, its 
coefficients are bounded by the maximum value of the 
function to be approximated. Besides, the magnitudes of 
the coefficients decrease rapidly. In General, the 
Chebyshev series has faster convergence rate than other 
series using orthogonal functions [12]. 

Assumptions 

A.l The time d u delay is known. 

A.2 A. ( q~ l ) and B n ( q~ l ) are coprime. 



3. Parameter Estimation 

The system described by equation (12) can also be 
written as: 

y t (k) = dj (j) i (& - 1) + cr (k) (13) 

where: 

4 (*-!) = (k ~ 1 )-y t (k - 2 ),...-y(k - n a ), 



P i {k) = P i {k-\)~ 

< 16 > 

4 (*) + 4 T (*)/>(*- lM(fc-i) 

with P { (0) = a 1 ,0 < a <oc and I is the identity 
matrix. 

A i {k) = A i0 A i {k-\) + (\-A i0 ) 
with A i0 = 0.95 and A, (0) = 0.96 . 



Mk) = 



y,(k) 

q.ik-l) 



Mk~ D = 



Mk-1) 

Vi ( k - 1) 



(17) 

(18) 



A.3 There exists a small positive scalar CO i such that 



\&i(k ) | 
V, (k) 



< CD i with q t is a normalizing signal. 



The proprieties which determine the parameter adaptation 
algorithm in adaptive control context are described in 
detail in [15]. 



4. Structure of the controller 



Together with the structure of Figure 1, the proposed 
control structure is given as: 



R t (q~ x )v ; (k) + G ( (q ~ l ) y. (k) = 

(19) 

A iJq-')yJk) 



where: 

v i {k) = u i (k) + u ic (k) et A im (q~ l ) is a stable 



m, (k - d u ), m, (k - d u - 1 (k - d u - n b ), 

T 0 (k),T l (k),...,T mA (k)] 

and 

6 t = [a iV a i2 ,...,a in b: o ,b* a ,...,b* inb , 



polynomial. 

U, (k) is synthesised from a PI control law using fuzzy 

logic technique and U ic ( k ) is the auxiliary control given 

in the sense to compensate the nonlinear effect, the cross- 
coupling and disturbances. 



h h h ] 

m is the expansion order of the Chebyshev polynomials. 
Based on parameterization (13), the identification 

algorithm giving estimates 6 t ( k ) of can be obtained 

using the normalized least-squares algorithm [13]. 

d i {k) = e i {k-\) + 



P,(k- lM(fc-l)g,-(fc)_ 

AW + tfik-DPtik-Dbik) 

e i (k) = y i (k)-0 i T (k-l)0 i (k-l) 



(14) 



(15) 




Figure 1 . Diagram schematic of the control 
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The problem includes the choice of the input and output 
variables of the fuzzy control rules, the generation of the 
rules and placing the membership functions. 

The fuzzy logic control law u t (k) considered in this 
paper is the form: 

Au i (k) = /(e^k^Ae^k)), i = l,...,n (20) 

Where <?, (k) = (k) - y im (k) , 

Ae i ( k ) = e i ( k ) - e i ( k - 1) and f: Rx R—> R is a 

mapping that realises the control strategy desired by 
membership functions and rules of the controller. 

Each input of the fuzzy logic controller (FLC) has tree 
fuzzy subsets: negative (N), zero (Z) and positive (P). 



Ri (q~ l )u, (k) + G, ( q~ l )y, (k) + 
R i (q- l )CjT(k) = A im (q- l )y im (k) 

with 

T du hjB^ 1) 

V * / -1 \ 



(24) 



(25) 



The polynomial R i (q 1 ) and G, (q 1 ) are given by: 



R i (q- 1 )A i (q- 1 ) + q- dii G i (q- 1 )B;(q- l ) = 

(26) 

A m (q~ l ) 




Figure 2. Membership functions for fuzzy sets of errors and change 
errors 



Where X- = {e t Aq ~ } and l x = |/ e / Ae |, parameters 

l e and / Ae define the widths of the membership 
functions as shown in figure 2. 

The output of the FLC is the change of the control signal 
A u i ( k ) . The final control signal is calculated by: 

u { ( k ) = u i (k- 1) + Au . ( k ) (21) 

For each FLC, the form of I th rule as the form: 

if e x ( k ) is N and Ae l ( k ) is P then A u\ ( k ) is A 1 

With A = [- Pvb >"Pb ’"Pm ’“Ps Aps »Pm >Pb »Pvb ]> 
where S signifies small, M medium, B big and VB very 
big and each element of the matrix A is the place of the 
singleton witch is used as membership functions. 
Weighted sum is used as a defuzzyfication method witch 
can be presented in a general form as: 

K 

A u t (k) = e^Pj/Uj (22) 

7=1 

Where jd ■ is the value of the membership function of the 
j th fuzzy subset of Au • (k) , K is the number of the 
fuzzy subset of Au i {k) , P. is the place of the j th 

subset in universe of discourse of Au • ( k ) and S is the 
scaling factor or fine tuning parameter. 

The control law V- ( k ) is given by : 



v i (k) = u i (k) + u ic (k) 

= u i (k) + CjT(k) 

Equation (19) can be written as: 



(23) 



The closed loop equation is given as follows: 

Operating on (12) by R i (q~ l ) gives: 

RM~ l ) A i (q~ 1 )y i (k) = 

q~ d "Ri (■ q~ X )B*i (q-')u, (k) + (27) 

Ri (q 1 T(k) + R[ (q~ l )cr (k ) 

Combining equations (24) and (27), we obtain: 

[Ri (q- 1 Mi 0T 1 ) + q~ du G, ( q - 1 )B: (q~ l )]y, (k) 

+ R i (q-')[q- d "Bl(q-')C i -^(t) = (28) 

q- d “ Bl (q-' )A im (q~ l )y im (k) + R t ( q~ l )*, (k) 
Combining equations (12), (24), (25) and (26) gives: 

4 m (q- { )[y, ( k ) - y im m = Ri (- q - l )°i (k) (29) 

It is easily proved in [13] that the boundedness of the 
input-output signal will follow from a sufficient 

smallness of the ratio < 0 ) i . It is noticing that, 

Vi (k) 

this assumption is less restrictive than the standard one 

H.(k) 

imposed on equation (7): is sufficiently small in 

V,(k) 

the mean. [13, 14]. 

Theorem : For the system (1), if the control law (24) is 
used under assumptions A. 1-3, then the closed-loop 
system is stable in the sense that: 

1. The output and all inputs are bounded for all the time. 

2. The output tracking error: ( k ) — y im ( k )} = 0 

k-> oo 

if <T (k) = 0 . 

The proof may be carried out along the same lines as in 
[16] and is then omitted. 



5. Simulation results 

Consider the two-input two-output (TITO) system having 
a polynomial nonlinearity described as follows [18]: 



A(q- l )y(k) = B d (q- l )Z(k) + Z(k) 



Where q 1 is the back shift operator. 
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y(k)eR 2 , Z(k) e R 2 and %(k) e R 2 are 
respectively the output, nonlinear input and the 
disturbances vector of the TITO system. The nonlinearity 
is assumed to be represented as follows: 

Z , (*) = fiO + fil U i (k) + - + fi P u r' (k) 

The TITO system can be decomposed into two MISO 
nonlinear subsystems as follows: 

Mq- 1 )y 1 (k) = q- d "B n (q- l )Z 1 (k) + 

q^B 12 (q- l )Z 2 (k) + ^(k) 

A ( q~ 1 )y 2 ( k ) = q~ dl 2 B 22 (q 1 )Z 2 (k) + 



q d2 ' B 2l (q 1 )Z l (k) + ^ 2 (k) 

With 

,, i n+0.35^ 1 +0.15o^ 2 ’ 

A(q "*)= . , 

1 + 0.72^ 1 +0.05g^ 2 _ 

, , r l + 0.85g _1 2 + 1.25g _1 

B, (q ) = y y 

Ll-12 + 0.65^- 1 1.65 + 0.23^ _1 

d n = 1 , d l2 = 2 , d 2l = 2 and d 22 = 1 . 

Zj (k) = 0.5«[ (k) + 0.25uf (k) 

Z 2 ( k ) = u 2 ( k ) + 0.83 w 2 ( k ) 

£ (k) = 0.1 rand (1 .600) 



14.35 -3.42 -3.20 10.5 0 

-5.3 10.47 -1.35 -14.35 
14.35 -3.42 -0.20 5.7 0 

-1.39 3.42 -0.33 -14.35] 



The number of rules for this example is 9 rules. 

It can be seen from Figures 3 and 4 that the performance 
of the decoupling via orthogonal polynomials and the 
proposed FLC is demonstrated and the tracking is 
obtained after small transient time. 

The fine tuning parameters was increased while the 
interactions were small enough. 
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Figure 3. The response of Ti (^) an< ^ 
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Figure 4. The response of y 2 {k ) and y 2m (&) 



g 2 (k) = 0.12 rand (1.600) 

The inputs are u l (k) and u 2 (k) . The outputs are 
y x (k) and y 2 {k) with % x {k) and % 2 (k) representing 
disturbances. 

After some manipulation, we obtain: 

A (A 1 )*(*) = q~ d "B* n (. q - 1 )«, (*) + //, (k) 

with q d " B u (q 1 ) = 0.5 q 1 +0.425^ 2 

A (A 1 )y 2 (*) = q dl1 b *22 (A 1 )« 2 W + H 2 (k) 

with q~ d21 B* 22 {q~ l ) = \.65q~ x + 0.25^ 2 



H l (k) and H 2 (k) are the function contains the 
nonlinear terms, cross-coupling and disturbances. 

To synthesise the FLC, the widths of the membership 

functions are: L = [l e , 1^,1^ , l ^ ] = [1, 1, 6.4, 6] . 

Let £ — 0.1 . 



The final values of L and £ 
simulation tests. 

The places of the membership 
consequence is chosen as: 



are obtained a few 
functions for all the 



The corresponding control laws are given on figures 5 
and 6. 
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Figure 5. The behaviour of the control signal (&) 
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Figure 6. The behaviour of the control signal U 2 (k ) 



6. Conclusion 

The approach presents a solution to the problem of robust 
control of MIMO non linear systems. For each MISO 
subsystem a local fuzzy PI control is given. The effect of 
nonlinear dynamics, disturbances and cross-coupling is 
compensated by projecting it on Chebyshev polynomials. 
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This leads to both solving the nonlinearity problem and 
to releasing the decoupling through the control synthesis. 
The proposed control allows the direct implementation 
for a large class of MIMO systems. The simulation 
results show the effectiveness of the fuzzy adaptive PI 
controller. 
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Abstract 

In this paper, a monitoring approach for fault detection in 
power system drives is presented. Due to its critical 
position in the power train the diagnosis of the two level 
three phase voltage inverter is of a great concern. Fault 
scenarios with single open-switch are considered because 
they are the most likely to occur. Several signals are 
analysed simultaneously in order to perform the 
diagnosis. The fault occurrence is revealed by a change in 
the mean value of a subset of the analysed signals. The 
diagnosis is realised in three steps. In the first step, 
signals are filtered using the SWT performed with the 
DB4 wavelet to extract the detail and approximation 
coefficients up to level 6. In the second step the 
approximation at level 6 is examined to detect changes in 
the mean. This is achieved with statistical hypothesis 
techniques. A Neyman-Pearson change in the mean 
detection test is used. Then at the third step, a signature 
table allows to isolate the faulty switch. The whole 
diagnostic procedure can perform on line because of its 
low computational cost. Real data recorded from a 
benchmark feed the proposed diagnostic tool. The results 
Presented here confirm the effectiveness of the proposed 
methodology. 

Keywords l Fault Detection and Isolation (FDI), 
Stationary Wavelet Transform (SWT), Neyman-Pearson 
detection test, power converter, hypothesis testing. 

1. Introduction 

Fault Detection and Isolation (FDI) is an essential aspect 
of process engineering. It is important not only from a 
safety viewpoint but also for the maintenance politics. 
FDI has received a considerable attention from industry 
and research communities because of the economic and 
safety impacts involved. There is however a number of 
practical challenges in designing robust FDI systems. 
Among the difficulties are the complexity of process, the 
lack of adequate models, the incomplete and uncertain 



data, the partial knowledge on the system, the effort and 
expertise required to develop and maintain the FDI 
systems. 

Signal processing is a well known tool to achieve FDI. It 
is used to analyze directly the signals measured online, 
avoiding system modelling. Its main drawback is that a 
change in some signal feature, e.g. the signal mean or the 
frequency contents, must be distinctive for each fault that 
may occur on the system. 

Detection tests that aim at detecting a change in the mean 
or the standard deviation of a signal are now very 
common [1, 2] techniques to perform diagnosis. 

Frequency representations are particularly useful for 
anomaly detection as for example in induction machines 
and drive systems. Indeed, extra frequency contents may 
appear under the influence of a particular fault. For 
instance, [3] deeply studies faults in a three-phase 
induction machine. The spectral analysis of electric and 
electromagnetic signals shows that mechanical 
abnormalities such as broken rotor bars generate 
characteristic frequency contents in the signals. 
Unfortunately, the Fourier Transform is unable to 
accurately analyse and represent a signal with non 
periodic features, such as transient ones. To study non 
stationary signals, time-frequency methods replace 
traditional spectral analysis [4, 5]. The Short Time 
Fourier Transform interpretation is close to a local Fast 
Fourier Transform analysis. The signal to analyse is 
multiplied by a sliding window with finite duration. Thus 
the spectrum is computed in real time and the variation of 
its contents is used to detect faults. 

The main drawback of STFT method is due to the 
constant time-frequency resolution, according to the 
Heisenberg-Gabor uncertainty principle. Indeed, there is 
a trade-off between time and frequency resolutions 
because an accurate time resolution requires a “short” 
analysis window while an accurate frequency resolution 
involves a “long” analysis window. However, in the 
context of diagnosis, this long analysis window 
introduces extra detection delays. 
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In order to obtain variable time and frequency resolutions 
with constant product, the Wavelet Transform (WT) has 
been introduced. Moreover, Wavelet analysis does not 
require stationarity hypothesis and it is well adapted to 
the analysis of signals with changes over the time. 
Multiresolution analysis (MRA) has been investigated for 
monitoring and diagnosis in various industrial areas [6, 7]. 
It is particularly suitable when no model of the physical 
phenomenon is available. In power converters, 
commutation failures and single phase short circuits at 
the AC side have been studied with such a tool. It has 
been shown that these faults produce time varying 
transients. [8] reports the use of the MRA to exhibit 
faults in power system drive. 

The Stationary Wavelet Transform (SWT) may be 
preferred to the Discrete Wavelet Transform (DWT) 
since it does not include any downsampling that may 
introduce a delay in the fault detection. The SWT fills in 
the gaps between the coefficients in any particular 
decimated DWT. 

In order to perform the FDI, the deviation of features 
from normal situation must be detected. Statistical 
hypothesis techniques can be applied in order to exhibit 
changes. From a practical viewpoint, the FDI framework 
usually implements two hypotheses, namely “unfaulty 
situation, or at least undetectable deviation from the 
normal behaviour” and “deviation from the normal 
behaviour”. It was shown in [9] that an optimal alarm 
system is fundamentally based upon a likelihood ratio 
criterion via the Neyman-Pearson lemma. The benefit 
from the trade-off between false alarm probability and 
missed detection probability is important to make 
decision rules. This allows to design an optimal alarm 
system that will elicit the smallest miss detection rate for 
a fixed false alarm rate. 

In this paper, the diagnosis of a power converter 
associated with an induction motor is considered because 
of its vital position in the power train. Faults in the 
converter may induce dramatic damages on the power 
system. Therefore, the development of a real-time 
diagnostic algorithm is of great interest. In order to study 
various fault scenarios, a simulator has been developed 
considering one single open-switch because these faults 
are the most likely to occur. Several signals must be 
analysed simultaneously in order to perform the fault 
isolation. 

The remainder of the paper is organized as follows: 
Section (2) focuses on the SWT used as a preprocessing 
tool to feed the detection and isolation phases of the FDI 
procedure. Section (3) emphasizes on the statistical 
hypothesis test that will be applied on the approximation 
coefficients during the detection step of the diagnostic 
procedure. Section (4) gives a brief description of the 
benchmark with its simulator and presents the diagnostic 
tool steps. The results achieved with the proposed 
method on real data are discussed in section (5). 

2. Stationary Wavelet Transform Theory 

A detection method based upon the MRA has been 
previously proposed in the context of power inverter 
diagnosis [10]. However, the main drawback of this non- 
redundant transform [14] is its non-invariance in 
time/space, i.e. the coefficients of a delayed signal are not 



a time-shifted version of those of the original signal. It is 
well admitted that time invariance is very important for 
diagnostic purpose and the DWT algorithm does not meet 
this requirement [15]. The SWT was introduced to make 
the wavelet decomposition time invariant. Among its 
appealing properties, the SWT improves the discontinuity 
detection capabilities of the DWT, as stated in [18]. 
However, the computational cost of the SWT compared 
with its MRA counterpart is significant, respectively O (N 
logAO and O(A0, where N is the length of the signal. Thus, 
it must be taken into account when online 
implementation of the FDI procedure is considered [17]. 
In the present paper, the SWT is used instead of the MRA 
for the detection of changes in the “low” frequency part 
of the signal while classical uses of the SWT in the 
diagnostic framework deal with the “high” frequency part 
of the signal [16]. The SWT method is used here for its 
filtering capability. It can be described as follows: at each 
level, when the lowpass and highpass filters are applied 
to the data, the two new sequences have the same length 
as the original sequence. To do this, the original data is 
not decimated. However, the filters at each level are 
modified by padding them out with zeros [11]. Suppose 
that a function f(x ) is projected at each step j on the 

subset V.(... cz V 3 cz V 2 c= cz V 0 ) . This projection is 

defined by the scalar product dj k of with the scaling 
function <fi(x) which is dilated and translated. 
a i,k = (l) 

0 m (jc) = 2-'#2-'jc-*) (2) 

where </>(x) is the scaling function, which is a low-pass 
filter, dj k is also called a discrete approximation at the 

resolution 2 j . If cp(x) is the wavelet function, the 
wavelet coefficients are obtained by 
d. k =(f(x),2- j (p(2- j x-k)) (3) 

dj k is called the discrete detail signal at the 

resolution 2 ] . As the scaling function (/>{x) has the 
following property: 

= Y h ( n )</>( x - «) 

4 4 n 

a j+l k can be obtained by direct computation from d f 
a j + l,k=Y h( - n ~ 2k ^ a i.n ( 4 ) 

n 

Moreover 

\<p(~b = Y s(n)(j>{x- n) 

4 4 n 

The scalar products ( f(x),2~ u+1) (p(2~ 0+1) x-k )) are 
computed with 

d j + l,k =Ys( n ~ 2k ) a i,n ( 5 ) 

n 

Equations (4) and (5) are the multiresolution algorithm of 
the traditional DWT. In this transform, a downsampling 
step is used to perform the transformation. That is, one 
point out of two is kept during the computation at one 
particular level. Therefore, the whole length of the 
function f(x) will reduce by half at each level. 
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However, for stationary or redundant transform, instead 
of downsampling, an upsampling procedure is carried out 
before performing filter convolution at each scale. Thus, 
the approximation coefficient a j+xk is obtained by: 

<W = 'L h i (n - k ') a i,n ( 6 ) 

n 

and the discrete wavelet coefficient is 

^j + u=TSj(n-k)a. n ( 7 ) 

n 

where hj and gj are the h and g filters with 2y-l-l zeros 
between each filter coefficient. The redundancy of this 
transform facilitates the identification of salient features 
in a signal. 

3. Neyman Pearson Approach 

When FDI is performed using signal processing 
techniques, the fault occurrence is supposed to change 
some of the signal characteristics. This change is 
exhibited with a fault indicator (the symptom) which can 
be built from a decision test established according to, for 
instance, statistical hypothesis tests. 

In binary hypothesis testing problem [1], the two 
hypotheses are denoted by H 0 and H { . The hypothesis H 0 
is known as the null hypothesis. This is because it usually 
represents the absence of fault (or at least undetectable 
fault). The hypothesis H x is known as the alternative 
hypothesis and it usually represents the presence of fault. 
P i is the probability distribution function of the signal 
under hypothesis Hi and P 0 is the probability distribution 
function of the signal under hypothesis H 0 . Two types of 
errors are defined. The first kind of error occurs when H 0 
is true and Hi is decided. This error is the “Type I error”. 
The probability corresponding to this type of error is 
known as the false alarm probability or false alarm rate 
a . The second kind of error occurs when H 0 is decided 
while Hi is true. This error is the “Type II error”. The 
probability corresponding to Type II error is named the 
miss detection probability (3 . 

The design of a hypothesis test involves a trade-off 
between a and (3 . The goal of detection theory is to 
choose an appropriate test or decision rule that is optimal 
in some sense: either it should minimize the probability 
of error, or it should maximize the probability of 
detection, or it should minimize the conditional risk, etc. 
Therefore, depending on the goal, there are different 
ways of customising the decision rule. The four situations 
are summarised in table 1 where: 
cr = Pr (H l I H 0 ) and f3 = Pr (H 0 1 H x ) 





Hq true 


H\ true 


Decided Hq 


D 00 =l-a 


£ 

ll 


Decided H \ 


D l0 = a 


D n =\-p 



Table 1. Retained decision 



In the literature, various detection tests can be found such 
as Bayes, Minimax, Neyman-Pearson... In the present 
work, the Neyman-Pearson test has been applied. This 
test maximises the detection probability (1-/3) if one 



accepts a certain probability of false alarm a . Fig. 1 
shows the probability distributions P () and Pi and the 



associated decision that corresponds to the threshold 
value A which is deduced from a . It also shows the false 
alarm probability and t a he miss detection probability (3 . 
It clearly exhibits the trade-off between a and [3 . 




The following test is now considered: 

Hr, : r = jUr, +1/ 

0 0 ( 8 ) 

H x \r = fii+v 

Therefore (8) corresponds to a change in the mean 
detection test, with ju 0 and ju x the means under 
hypotheses H 0 and Hi respectively, v is a white noise. 
Let P x = Pr(rl//,) and P 0 = Pr(rl/7 0 ) . The Neyman- 
Pearson lemma involves the so-called likelihood ratio 
V(r) defined by [1]: 



T „ , Pr mo ? 

V(r) = — < rj 

Pr (r\H 0 ) 

If V(r) passes a threshold f], the choice is H u otherwise 
the choice is H 0 . The threshold is computed under the 
constraint on the false alarm probability: 



a 



= J Pr(r,ltf 0 )dr 



The usual way is to consider the logarithm of the 
likelihood ratio, named the log-likelihood ratio: 



W(r) = ln(V(r))> In (77) 



(9) 



Under the Gaussian noise hypothesis ( v is N(0,a ) ), it 

can be shown that: 

2 



A = - 



M\ -Mo 



-111(77)+ 



M 1 -M 0 



( 10 ) 



In this work, the case of a multidimensional signal is 
considered. The hypotheses become: 

H 0- r n=Mo+ V 'n’ n = l:N 

H l :r n =ju l+ v'„; n = l:N 

The Gaussian white noise independent identically 
distributed (iid) case is obvious and it is not presented in 
this paper. From now on, the noise sequence { v n '} is 
supposed correlated with covariance matrix C . This 
situation may arise in linear filtering, r n being computed 



with: 
r n =^X 



( 12 ) 
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where h = [/Zp---,/z L ] r e $R L is the filter impulse response 

and x n = [ W * • • , x n _ 1 ,x n f e . Consider that the raw 
data x n satisfies: 

x„=x n +v„ (13) 

x n is the deterministic part of the sample x n . v n is 
supposed to be a zero mean Gaussian white noise 
A(0, cr 2 ) and the sequence {v n } are iid variables. 

Therefore: 
r n =h r x„+h r v„ 

= r n +y' n 

and 

x,=[W"*J eSR 1 . (15) 

The elements of v \ = [v \_ L+l • • • v \ f e 91 l are correlated 
noises with covariance matrix equal to: 

C = a 2 AA T (16) 

where 



(14) 



A = 






K o 



o 



0 }\ h 2 ••• h L 0 0 

0 0 0 \ h 2 ••• h L 



: ^Vx(A^+L-l) (17) 



Without lack of generality, consider ju 0 =0 in (11). The 
hypotheses become: 

H 0 :r n =v' n - n = l:N 

H i ■ r n=M l +y'n', n = \:N 

Define r = \r v ---,r N ^ e 9{ N and m, = [//,, — e 5R W . 

Then, the probability functions under hypotheses H 0 and 
Hi are respectively: 



P 0 =Pr(rltf 0 ): 
P x - Pr(rl// ] ) = 



1 



<J(2 7r) N det(C) 

1 y (( r_m i f c 1 ( r_m i )) 



^(2 7r) N det(C) 
The likelihood ratio becomes: 

1 



(18) 

(19) 



V(r) = 



-\l(27r) N det(C) 



y(( r_ m 1 ) T C 1 (r— mj )) 



V(r) = e 



yj(2 7r) N det(C) 

— l^-2r T C _1 m 1 +mic _1 m 1 ) 



y((rfc-‘(r)) 



( 20 ) 



and its logarithm (9) is: 

ln(V(r)) = — (-2r r c _1 m, +mfc _1 m 1 ) ^ ln(77) 

2 V Hq 

«r r C m ( > ln(t7)H — — L 

Z 

As expected, the noise correlation must be taken into 
account and obviously, the threshold depends on matrix 
C. Note that (20) corresponds to the computation of a 



weighted sum of the samples r while the iid case 
considers the unweighted sum ^ ^ r . 



Even if the computation of (20) is of no difficulty, a 
suboptimal solution is derived with the noise sequence 
{v' n } supposed to be uncorrelated. Consider the signal 
mean: 



Y n 

m = — V r 

N i n 

1 v n= 1 



Thus, the hypotheses to be tested are: 



H 0 \ m- v" 



( 21 ) 



( 22 ) 



H x : m = + v " 

where v" is supposed to be a Gaussian white noise 
N(0,<r' a ) and 



cr = — 
N 



?{[' ' - T ! 



(23) 



The likelihood ratio is given by: 



V(m) = 



V(2 ~n)o" 



1 

=e la 

.2 






V ( 2;r ) cr " 2 



<^> V (m) = e 






7 



It follows that the detection test to be implemented is: 

Ai 



a-" 2 



m \ 2 = 

HQ 



\n{r 1 '\ + — 

A V ^ 2 



(24) 



This threshold is similar to the one computed in (10). 
This suboptimal detection test is implemented in the FDI 
procedure, the mean value of the deterministic part of the 
analysed signal over the window of length N being null. 

4. Benchmark description 

As stated above, the power converter reliability is 
extremely important in industrial applications. Thus a 
large amount of papers deal with the faulty mode 
behaviour of static converters, and with the fault tolerant 
control and protection of voltage source inverter systems. 
Some of them consider various fault modes of a voltage 
source PWM inverter for induction motor drive [12]. 
Others are interested in fault tolerant control of induction 
motor drive applications using analytical redundancy, 
providing solutions to the most frequent faults [13]. 

This paper deals with the fault detection and isolation 
problem of a voltage source PWM motor drive system. 
The one open switch fault case is considered, due to its 
frequent occurrence. A benchmark depicted in Fig. 2 has 
been developed. It consists of a three-phase inverter with 
IGBTs and anti-parallel diodes controlled by a DSpace 
DS1103 board inserted into a PC Pentium and an 
induction machine as a load. This mixed RISC/DSP 
digital controller, based on two microprocessors (power 
PC 604e-333MHz and TMS320F240-20MHz) and four 
high-resolution analog-to-digital (A/D) converters (0.8 ps 
-12 bits), provides an efficient platform for real-time 
floating point calculation. The controller board is fully 
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supported by the Matlab/Simulink environment. The 
main voltage supply is set to 200 V. The characteristics 
of the induction motor connected to the inverter are 1.8 
kW, 380V, 4.5A and 1460 rpm. 

The electrical drive has been tested under normal 
conditions and with an open circuit fault on the IGBT T4. 
Fig. 3 shows data recorded on the setup. The three phase 
stator currents and the control of IGBT T4 waveforms are 
shown. After the fault occurrence, the phase 2 current has 
a positive DC offset and a distorted waveform. This is 
due to the fact that during a commutation in leg 2, the 
load in phase 2 is connected to the DC positive voltage 
bus all the time through the defective IGBT. Negative 
DC offsets are then present in both phase 1 and 2 currents. 



/Inverter with integrated sensors^ ^ 




(b) 

Fig. 2 (a) Configuration of the experimental setup, (b) experimental 
setup 



As previously mentioned, faults are common phenomena 
in power inverter. They induce device destruction. Due to 
their cost constraints, a power inverter simulator using 
IGBT and anti-parallel diode detailed models was 
developed in the Matlab environment. The simulator has 
been set in the same operating conditions as the 
laboratory setup. Parameters are set equal to those of the 
benchmark. Fig. 4 exhibits the results obtained with the 
Matlab simulator. Simulation results have been compared 
successfully with measurements in steady state and 
transient conditions. It can be noticed the excellent 
agreement between the simulated and the recorded 
signals in both normal and faulty conditions [8]. 




10 

~ 0 
-10 

1.55 1.6 1.65 1.7 1.75 

t(s) 

Fig. 4 Simulation results, zoom around the fault time occurrence, (a) 
control of IGBT T4, (b) Phase 2 stator current, (c) Phase 3 stator current 
and (d) Phase 1 stator current 

5. Fault Detection and Isolation 

The current signals in the three phases are analysed 
simultaneously. The output of the detection phase 
provides six symptoms that feed the isolation step. The 
fault isolation is performed with a signature table to 
isolate the faulty switch. Therefore the FDI procedure is 
organised in three steps Fig. 5: 

- Preprocessing of the stator currents with the SWT to 
obtain the approximation a 6 at level 6; 

- detection using Neyman Pearson tests (change in the 
mean detection tests); 

- isolation with a signature table. 





Chi f 100 V |Ch2| 2.00 V |M |20.0ms| A| Ch2 X -1.08 V 

SSEI 2.00 V ^Ch4| 2.00 V 



Fig. 3 Switch-on failure transient waveform: (1) control of IGBT T4, (2) 
phase_2 stator current, (4) Phase_3 stator current and (3) phase_l stator 



current 




Fig. 5 FDI blocks diagram 



In this section the three steps are presented. Then the 
results obtained with the data recorded from the 
benchmark are discussed. 



5.1 SWT Preprocessing 

The SWT is applied to the three- phase current signals 
using the Daubechies wavelet DB4. It decomposes the 
signal into approximation and detail coefficients. The 
signal analysis was performed up to the sixth level. 
Indeed, the approximation at level 6 emphasises the fault 
occurrence while high frequency contents are sufficiently 
filtered. A single transistor base drive open circuit of the 
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power converter is considered. An extensive 
investigation has been carried out for various faulty 
IGBTs. The mean value (on a time period) of the stator 
currents in steady state is zero. When the open circuit 
fault appears in the inverter transistor, the fault 
information affects several frequency bands Fig. 6. 



D1 




t(s) 



Fig. 6 Six levels of detail and approximation of phase 1 stator currents 
before and after open-circuit fault on IGBT T1 

Spikes in the details at the sixth level and in the 
approximations can be seen. Their appearance time is 
directly correlated with the occurrence time of the fault. 
Moreover, this kind of fault introduces a non-zero DC 
offset that can be clearly observed in the approximation 
at level six. Fig 7 depicts a current and its approximation 
a 6 at level 6 when IGBT T1 fails. Fig. 8 shows that a 
larger DC offset appears in the faulty phase when 
compared to the other phases. 

As a result of time domain studies, the faulty situation 
can be discriminated from the healthy one. Note that the 
6th level approximation a 6 is greater than its DWT 
counterpart. Therefore, the gap (i.e. the mean change 
between the healthy and the faulty situations) is more 
easily detected in an automatic way. As a consequence, 
the detection of the mean change in a 6 is facilitated 
compared to the detection of the change in the current. 
However, a small unbalance of the load must not be 
confused with the faulty scenarios. Thus, the threshold 
has been computed in order not to introduce false alarms 
due to unbalance of the load. 

5.2 Fault detection with Neyman Pearson 

Approximation at level 6 is analysed with the Neyman- 
Pearson test presented in section 3. The false alarm 
probability a is set to 0.05 and the miss alarm probability 
is minimized. 

The measured currents are sampled at T e = 0.2 ms. These 

raw data are supposed to be modelled with: 

i(t) = i(t) + v(t), t = kT (25) 

where i is a periodic deterministic signal, and v is a 
Gaussian white noise A(0,cr 2 ) . When a fault occurs, the 
current is supposed to be: 

i(t) - i ( t ) + fi x + v(t) (26) 

with ju x a constant. 

The SWT is applied to i and the approximation a 6 is 
considered. Actually, the SWT is a linear transformation. 



Therefore, with respect to (11), in normal condition, one 
gets: 



aAt) = W\ = aAt) + v'(t) 



(27) 



with a 6 a periodic deterministic signal and v ' a 
Gaussian noise with zero mean and cr 2 H T H variance. 
Note that the sequence { v ’} is no longer independent. Its 

covariance matrix is given by (16). 

The Neyman-Pearson test is performed on m (21): 

1 N 

mW = 77ZA(/-* +1 ) 

A /=! 

= 6 (t-i + 1) +— yV^-i' + l) 

The window length N is chosen so that 
= The choice of N- 200 samples 

corresponds to 2 signal periods. Therefore, the 
hypotheses (19) are considered: 

H 0 :m(t) = v"(t) 



H ] : m(t) = ju x + v"(t) 
where v" is defined in (20). 

For each stator current, two tests are implemented. 
Indeed, the sign of the mean change in a 6 is exploited to 
isolate the fault. The minimum mean change to be 
detected ju x is chosen so that the test is insensitive to a 



small unbalance of the load. The mean ju x is set to + 20 
or -20: the sign of the change will be used for isolation 
purpose. The test for the positive (negative) mean change 
is named T u (, a6 ij) (T d (a6 ij )), 7 = 1:3. The six Neyman- 

Pearson tests are then derived and their thresholds are 
computed in accordance with (21). For each test, the 
symptom has a binary value: “0” means that “hypothesis 
H 0 is kept” while “1” stands for “hypothesis H x is 
selected”. 

Figs. 7 to 10 depict two scenarios which correspond to 
the open switch on IGBT and T 4 respectively. In the 
experimentations the load is balanced. Thus, the currents 
have the same amplitude in normal operating mode. 
Moreover, a transient phase can be seen from time 0.2s 
until time 1.35s. 

Fig. 7 shows the stator current ii in phase 1 (blue curve) 
in the case of a fault on IGBT T i and its approximation at 
level 6 (red curve). 
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Fig. 7 Stator current in phase 1 (il) for the fault on IGBT T1 and its 6th 
approximation <z 6 
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Fig. 8 depicts the three phase stator currents for the 
related fault scenario (top of the figure) and the six 
symptoms computed with the Neyman-Pearson tests. As 
expected, the detection tests T d (a6 q), T u (a6 i 2 ) and T u (a6 
z 3 ) are sensitive to this fault. 

Fig. 9 shows the stator current i\ in phase 1 (blue curve) 
for the second fault scenario (faulty IGBT T4) and its 
approximation at level 6 (red curve). 

Fig. 10 depicts the three phase stator currents (top of the 
figure) and the six symptoms obtained from the Neyman- 
Pearson tests. For this scenario, the detection tests T d (a6 
z'i), T u (a6 i2 ) and T d (a6 i 3 ) are sensitive to this fault. Note 
that the mean changes are different for both scenarios. 







Time [s] 

Fig. 8 Three phase stator currents for the fault on IGBT Ti (top) and 
symptoms derived from the six Neyman-Pearson 
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Fig. 9 Stator current in phase 1 (il) for the fault on IGBT T 4 and its 6th 
approximation a ( . 
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Fig. 10 Three phase stator currents for the fault on IGBT T 4 (top) and 
symptoms derived from the six Neyman-Pearson 



5.3 Fault Isolation 

The signature table was obtained with the simulator by 
making an exhaustive study of the various fault scenarios. 
The test results presented in section 5.2 make possible the 
isolation of the faulty switch. Indeed each combination of 
symptoms characterizes a fault. Thus the signature table 
is strongly locating. Table 2 gives the signature table. “0” 
means that hypothesis H 0 has been chosen for the 
considered test while “1” stands for H x has been kept. 

The isolation table has been used for the diagnosis of the 
setup plant. The fault isolation is based on a set of rules 
that are deduced from the isolation table. 



Detection 

Test 


Fault 

on 

IGBT ! 


Fault 

on 

igbt 2 


Fault 

on 

igbt 3 


Fault 

on 

igbt 4 


Fault 

on 

IGBT S 


Fault 

on 

igbt 6 


T u (a6h) 


0 


1 


1 


0 


1 


0 


T d (a6 i,) 


1 


0 


0 


1 


0 


1 


T u (a6 i 2 ) 


1 


0 


0 


1 


1 


0 


T d (a6 i 2 ) 


0 


1 


1 


0 


0 


1 


T u (a6 i 3 ) 


1 


0 


1 


0 


0 


1 


T d (a6 i 3 ) 


0 


1 


0 


1 


1 


0 



Table 2. Signature table to isolate the faulty IGBT 

Two scenarios have been considered on the real plant, 
namely, a fault on Ti and a fault on T 4 . The results 
presented in Figs. 7 to 10 are coherent with the simulated 
ones. Therefore, it is concluded that the single open- 
switch fault for all the possible scenarios can be detected 
and that the faulty switch can be isolated. 

6. Conclusion 

This paper presents a fault detection and isolation 
technique based on three steps. The first one makes use 
of the Stationary Wavelet Transform (SWT) as a 
preprocessing tool. The second step implements a 
statistical hypothesis testing method to detect changes in 
the wavelet coefficients. A signature table is used to 
isolate the fault during the third step. 

The SWT has been chosen as an alternative tool to the 
Discrete Wavelet Transform (DWT) mainly for its time- 
invariant property and its filtering capabilities. The 
approximation at a given level is analysed in the second 
step of the procedure. A Neyman-Pearson detection test 
is implemented to detect the fault occurrence. This test 
takes into account the statistical characteristics of the 
filtered noise. Indeed, this latter cannot be any more 
considered as an iid sequence. Several detection tests run 
in parallel, providing symptoms that are used with the 
signature table to isolate the fault. 

The proposed FDI technique has been successfully 
applied to a power inverter associated with an induction 
machine which is widely spread in industrial applications 
and a basic part of electric vehicles. The power inverter is 
a critical element of the power train: a fault on a switch 
may induce serious damages on both the motor and the 
inverter itself. The single open switch fault case was 
considered for the various IGBTs. The three phase stator 
currents are simultaneously analysed with the SWT up to 
level 6. Each one of the three approximation coefficients 
is processed in parallel by two Neyman-Pearson 
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detection test that compute symptoms. Then the six 
symptoms allow the isolation of the faulty switch. 

A simulator has been developed in order to carry out 
extensive investigations of the various fault scenarios 
which have been exploited to define the signature table. 
Real data recorded on the plant have been used to 
validate the FDI algorithm. 
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Abstract 

This paper deals with the modeling and computation of 
rotational hysteresis core losses in a rotating surface- 
mounted permanent magnet (SMPM) machine. 
Formulations of these losses taking into account the 
rotational behaviour of the magnetic field are developed. 
These formulations are used together with the time 
stepped finite element method to calculate the hysteresis 
core losses in the considered design according to the 
distribution of the rotating magnetic field. The originality 
of this work resides on the search of the optimum number 
of directions to be used for Vector Preisach model. This 
parameter has an interesting impact on the magnetic 
intensity value, and hence on the hysteresis loss value. So 
it seems significant to optimize this parameter by 
considering the two main criteria: (i) accurate estimation 
of the hysteresis iron losses, (ii) reasonable computation 
time. 

In a first step, the effect of the number of directions on 
the hysteresis loss term has been investigated. In a 
second step, the hysteresis losses, considering the 
optimum number of directions and the developed 
formulations have been investigated. The calculated 
results are discussed. 



a angle between B and H 

Iff I magnetic field intensity magnitude 

\ L± e,ti | 

|# | induction magnitude 

H exti x- of H 

H e y, tl y-of/7 

B ex>ti x- of B 
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e element 

N Total number of stepped FEA steps 

p Pair pole number 

L Stack length of the stator 
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B m i Peak induction of the i th element 
e Unit vector along the direction 6 n 

6 n Polar angle 

A 6 Discritization angle 
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R FP Fixed point residual 

N d Number of directions 



Keywords: magnetic hysteresis losses , rotating electrical 
machines, time-stepped finite element analysis. 
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Nomenclature 

frequency of the magnetic flux waveform 
peak of the sinusoidal flux density 
Steinmetz constant 
hysteresis constant 
constant 

angle between B and the x-axis 



1. Introduction 

A major part of the losses in electrical machines is the 
losses in the iron core. For an effective design of electric 
machines it is important to get accurate information about 
these losses. This is possible thanks to the knowledge of 
the magnetic field in the machine. 

Because of the complicated geometry and operation 
modes of electric machines, the magnetic field variation 
is complex; hence a numerical technique is needed to 
solve this matter. In this way, the finite element analysis 
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(FEA) has been used to compute the magnetic field and 
thus the iron losses. 

Generally, the iron losses are separated into hysteresis 
loss, classical eddy currents loss and excess loss. The 
present work deals with the computation of the first iron 
loss term. 

In this paper an investigation of the hysteresis iron losses 
at no-load operation in a surface-mounted permanent 
magnet motor has been outlined. 

Generally, in time- stepped FEA of electrical machines, 
single-valued curves are usually used for all materials 
comprised in the model. In this case, the irreversible 
material behavior in the iron cores completely ignored in 
the time stepped magnetic field equations, but the iron 
losses may readily be predicted a posteriori [1]— [3] . 
Another way to efficiently cure the negligence of the 
irreversible material behavior is the inclusion of an 
adequate hysteresis model in the FEA. In our case the 
used hysteresis model is the Vector Preisach Model, 
which has been incorporated in FEA in a reversed 
fashion. The vector Preisach model is defined as a 
superposition of scalar Preisach models regularly 
distributed along all possible angular directions [4] (a 
scalar model along each direction on the space). In 
practice, the number of direction is taken finite. This 
affects directly the value of the magnetic field H, and 
hence the hysteresis losses value. For these reasons, it has 
been the subject of deep investigations in the present 
work, in an attempt to reach the optimized number of 
directions to be adopted in the VPM used in this work. 
The considered design is a surface mounted permanent 
magnet motor (SMPM). In order to reduce the 
computation time, the calculation is made up in four 
regions of the considered topology. 

2. Hysteresis loss formulations 

According to the statistical theory [5], the average power 
loss per unit volume consists of the sum of the hysteresis, 
classical and excess loss contributions. In the present 
work, the hysteresis loss is the subject of investigation. 
Many empirical and analytical formulae, dealing with the 
hysteresis loss modeling, depend on the nature of the flux 
density variation have been proposed in the literature. In 
a rotating machine, owing to the slotting and saturation 
effects, the magnetic flux density may be significantly 
distorted from the ideal sinusoidal variation. In addition, 
in some parts of the machine the magnetic flux density 
becomes rotational rather than alternating. Hence, a 
hysteresis model should take into consideration this 
rotational behavior of the magnetic flux density. For a 
long time, the hysteresis losses in rotating electric 
machines have been computed using alternating core loss 
model, in which the magnetic flux density is supposed to 
be purely alternating. This is due essentially to the lack of 
the data associated with the rotational core loss and the 



lack of an appropriate model According to the Steinmetz 
equation, the measurement and the calculation of the 
hysteresis losses are normally carried out with sinusoidal 
flux density of varying magnitude and frequency. These 
measurements and calculations are based on the 
following formula 

P h= k hf B J 3 (D 

/ frequency of the magnetic flux waveform, 

B m peak of the sinusoidal flux density, 

/? Steinmetz constant 

k h hysteresis constant 

For correction and modification, various models have 
been proposed [6] -[8], 

p h =hfB,! a+bP) ( 2 ) 

Where a and b are constant. 

In a permanent magnet motor the magnetic flux density is 
far to be purely alternating, in certain region the flux 
density is purely rotational. So, a hysteresis loss model 
should take into account this rotational behavior of the 
flux density. For the modeling of rotational hysteresis 
loss, a basic Stoner-Wohlfarth model has been proposed 
in [9] and an empirical formula has been developed in 
[10]. Furthermore, the numerical implementation of such 
models is totally impractical. 

Pfi = F halt + Phrot = ^ C (Phalt ) (^) 

In [11] an estimation of the effect of rotational losses in 
an induction machine has been proposed. This approach 
estimates thee iron losses by means three method: 
considering the losses provided by: purely alternating 
field, due to the orthogonal components of the flux 
density and by applying a correction factor based on 
experimental data to improve the rotational loss 
calculation. In this case it has been shown that the third 
method provides more accurate results than the other 
two. The lack of this method is the introduction of a 
correction coefficient depending on the flux density peak; 
this would considerably increases the computation time. 
An investigation of the iron losses modeling has been 
reported in [12], in which the calculation of the iron loss 
in rotating electrical machines by using different iron loss 
models has been studied. The accuracy of each iron loss 
model was identified by the comparison between 
measured and calculated iron losses in the case of 
distorted elliptically flux density waveforms. In these 
models the calculation of the iron losses is based on the 
equivalent orthogonal alternating flux density 
components which. 

An interesting iron loss model has been proposed in [13, 
14]. In this model the hysteresis losses are expressed as 
follows: 
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Where T is the period of the magnetic field and p h (t) is 
the hysteresis power loss. 

According to [14] the hysteresis loss in (4) can be 
expressed as follows: 



Pk (0 = ^ J I# |^cos(a)fr ■ + ^ J ™{H x B\dt (5) 



T i dt 



Where 0 is the angle between the magnetic flux field B 
and the x-axis and a is the angle between B and H. The 
first term corresponds to the loss under purely alternating 
dd 

induction, in which is equal to zero and hence it can 

dt 

be expressed as follows: 

d\ B \ 

TV dt 



P halt (0 — T \ 1^1 



-dt 



( 6 ) 



Whereas the second term in (5) corresponds to the loss 

dB 

with purely alternating magnetic field ( — = 0 ) and it 

dt 



can be expressed as follows: 

p hr jt)=vV (HxB )^t 

T i dt 



(7) 



This model will be adopted in the present work. And the 
total hysteresis power loss is obtained by summing the 
alternating and rotational components: 



1 , ,d\B\ 

Ph = Pah + Phrot( t ) = —\\ H \ ~r~ c 
T t dt 



f + -j—(HxB) dt (8) 
T t dt V h 



By performing a 2-dimensional time-stepped FEA, the 
expression in (8) is transformed as follows: 
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Where p , L, A h N e , B mi are respectively the pair pole 
number, the stack length of the stator, the area of the 
finite e element i and the total number of finite element 
mesh. The last equation has been coupled with the finite 
element analysis to calculate the hysteresis iron losses in 
a surface mounted permanent magnet motor. Owing to 
the symmetry and periodicity in the structure, the 
problem has been reduced to only one pole of the 
considered topology (figure 1). Doing so, the value of the 
obtained hysteresis losses is 34 2 W. Which, when 
compared to the one obtained is [15] (same value) shows 
that such a model is satisfactory. 




Figure 1: Fourth of the studied topology [16] 




d\B. 



dt 



'-dt + , ^ 

1 t dt 



Where e t - 1 , ^B e t - 1 , H ex ,ti > H e x,ti> P ex, ti ctnd B e y t i are 

respectively the magnetic field intensity and induction 
magnitudes, and the orthogonal components respectively 

of H and B. f=l N, N is the total number of stepped 

FEA steps. From the FEA the flux density and the field 
intensity for each element can be determined, taking into 
account the symmetry and periodicity in the structure. 
The analysis is limited to a half period. Then, the interval 
between consecutive time steps is: 

A t = t i -t i _ l =— (10) 



The magnetic hysteresis losses per element are derived 
from equation (9), then by summing these losses over the 
elements number the total hysteresis losses can be 
obtained as follows [17]: 



3. Hysteresis Modelling 

In this paper, the hysteresis in ferromagnetic parts (stator 
core) is handled by means of the vector Preisach model, 
consisting of angularly distributed scalar models, which 
is applied in a reversed manner. 

Considering a finite number of directions N d equally 
spaced as illustrated in figure 2, H is computed as [18]: 

2 Nd 

tf(/) = -A#£ H n sv (B£ ai )e 6h ( 12 ) 

^ „= 1 

e fo is the unit vector along the direction specified by 

sv 

polar angle 0 n , AB is the discritization angle and H n 

is the scalar Preisach corresponding to the direction n. 

As presented in equation 12, the value of the magnetic 
field intensity H depends directly on the number of 
directions used in the vector Preisach model. This 
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influences explicitly the value of the hysteresis loss term 
given in equation 1 1 . 




Figure 2: Finite number of direction. 

For the integration of VPM in FEA an interface was 
needed. This interface is obtained by means of the fixed- 
point technique iterative procedure, in which the 
relationship between B and H is written as follows [17], 
[18]: 

B = juH + R fp (13) 

Where p is the fixed-point coefficient and R FP is a 
residual founded iteratively. Coupling equation (13) to 
appropriate Maxwell equations, the magnetic problem is 
formulated in terms of magnetic potential vector. The 
problem solution starts with an arbitrary given value of 
R fp , providing the magnetic potential vector A, from 
which B is calculated and supplied as input of VPM 
giving the value of H. With B and H known, a new value 
of R fp is calculated using formula (13). This value is 
again supplied in FEA equations. The procedure is 
repeated until the problem convergence is reached. 

4. Effect of the Number of Directions 

As discussed in the previous section, the number of 
directions used in the VPM has an increasing impact of 
the hysteresis losses term, within this intention a deep 
2D FEA based investigations have been carried out in 
order to illustrate the effect of this parameter (number of 
direction) on the hysteresis losses. Doing so, the 
hysteresis losses in different elements have been 
computed for different directions. The obtained results 
are illustrated in figure 3. 
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Figure 3: Hysteresis losses versus directions number 




Figure 4: Points for which the hysteresis losses has been 
computed 

As illustrated in figure 3 one can notice that the value of 
the hysteresis losses, for all regions, becomes slightly 
constant for a number of directions corresponding to 4. 
Then this parameter should be carefully chosen because 
low values lead to a bad estimation of the value of the 
iron losses and high values lead to the increase of the 
computation time. 

In addition to the number of the VPM directions, there 
are many other parameters affecting directly the value of 
the hysteresis losses, such as the machine geometry 
essentially the pole number. For this reason, the effect of 
this parameter on the hysteresis losses has been 
investigated; the obtained results are illustrated in figure 
5. 
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Pole number 




Figure 8: Flux lines plot: 6 poles 



Figure 5: Hysteresis losses versus pole number 

As illustrated in figure 5, one can clearly notice the 
increasing of the hysteresis losses with the pole number; 
this leads to the reduction of this parameter to obtain a 
low value of the hysteresis iron losses. However, the 
studied machine is a synchronous concept where the 
output power is proportional to the pole number, so a 
reduction of this parameter leads to the output power 
reduction. In this case a compromise should be taken. 
The flux lines plots corresponding to each pole number 
are illustrated in the follows figures. 

Referring to these figures, we can notice on one hand the 
increase of the flux leakage and on the other the increase 
the field variation with the pole number yielding to a 
hysteresis losses increasing. 

The effect of the other geometric parameter is implicitly 
considered in the stator teeth and yoke areas, witch 
appears in the hysteresis losses expression in formula 1 1 . 




Figure 6: Flux lines plot: 2 poles 




Figure 7 : Flux lines plot: 4 poles 





5. Conclusion 

The hysteresis iron loss investigation in a SMPM 
machine at rotational field operation has been presented. 
For this reason, an analytical model was developed, and 
it has been used in conjunction with time-stepped FEA to 
estimate the hysteresis iron losses. This model provides 
reasonable accuracy, because it takes into consideration 
the rotational behavior of the magnetic field. 

A brief study of the inclusion of the VPM in time stepped 
FEA has been also investigated. In this case a special 
attention is given to the effect of the number of directions 
used in the VPM in the hysteresis losses. The 
investigation of this parameter is required by the fact that 
in the expression of the hysteresis losses appears the 
magnetic field intensity H witch depends clearly on the 
direction number. In this work, it was shown that this 
parameter has a crucial impact on the hysteresis losses 
value: for low values of this parameter the hysteresis 
losses value varies. Whereas, for high values the 
hysteresis loss becomes almost constant, furthermore the 
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computation time would also increases and hence a 
compromise should be taken. In this case a number of 4 
directions were taken and would be adopted in order to 
accurately estimate the hysteresis losses. The effect of 
the pole number on the hysteresis losses has been 
investigated. In this case it has been shown the increasing 
of the hysteresis losses with this parameter; and hence the 
hysteresis losses reduction requires the reduction of this 
parameter. However, the studied machine is synchronous 
concept where the output power is proportional to the 
pole number, so a reduction of this parameter leads to the 
output power reduction. In this case a compromise should 
be taken 

References 

[1] G. Bertotti, A. Boglietti, M. Chiampi, D. 
Chiarabaglio, F. Fiorillo M. Lazzari, “An improved 
estimation of iron losses in rotating machines,” 
IEEE Transactions on Magnetics , Vol. 27, pp. 5007- 
5009, November 1991. 

[2] L. Dupre, R. Van Keer and J. Melkebeek “An iron 
loss model for electrical machine using Preisach 
theory,” IEEE Transactions on Magnetics , Vol. 33, 
pp. 4158-4160, Sept. 1997. 

[3] J. C. Gyselinck, L. Dupre, L. Vandelve and J. 
Melkebeek “Calculation of no load induction motor 
core losses using the rate-dependent Preisach 
model,” IEEE Transactions on Magnetics , Vol. 24, 
no. 4, pp. 621-630, Nov. 1998. 

[4] J. Saitz “Computation of the Core Loss in an 
Induction Motor Using the Vector Preisach 
Hysteresis Model Incorporated in Finite Element 
Analysis,” IEEE Transactions on Magnetics , Vol. 
36, no. 4, pp. 769-773, July. 2000. 

[5] G. Bertotti, “General properties of power loss in soft 
ferromagnetic materials,” IEEE Transactions on 
Magnetics , Vol. 27, pp. 5007-5009, 1989. 

[6] M. Amar and F.Protat, “A simple method for the 
estimation of power losses in silicon iron sheets 
under alternating pulse voltage excitation,” IEEE 
Trans. Magn. vol. 30, pp. 842-944, March 1994. 

[7] G. Bertotti, F. Fiorillo and M. Pasquale, 
“Measurement and prediction of dynamic loop 
shapes and power losses in soft magnetic materials,” 
IEEE Trans. Magn. Vol. 29, pp. 3496-3498, Nov. 
1993. 

[8] G. Bertotti, “Dynamic generation of the scalar 
Preisach model of hyteresis,” IEEE Trans. Magn. 
vol. 28, pp. 2599-2601, July 1992. 

[9] H. Pfutzner, P.schonhuber, B. Erbil, G.Harasko and 
T. Klinger, “Problems of loss separation for 



crystalline and consolidated amorphous soft 
magnetic materials,” IEEE Trans. Magn. vol. 27, pp. 
3426-3431, July 1991. 

[10] W. F. Archenhold, H. F. Sandham, and J. E. 
Thompson, “Rotational hysteresis loss in grain- 
oriented silicon-iron,” British J. Appl. Phys ., vol. 11, 
pp. 46-49, Jan. 1960. 

[11] R. D. Strattant and F. J. Young, “Iron losses in 
elliptical rotating fields,”/. Appl. Phys., vol. 33, no. 
3, pp. 1285-1286, Mar. 1962. 

[12] Carlos A. Hernandez- Aramburo, Tim C. Green, and 
Alexander C. “Estimating Rotational Iron Losses in 
an Induction Machine,” IEEE Trans. Magn. vol. 39, 
pp. 3527-3533, 2003. 

[13] B. Stumberger, V. Gorican, G. Stumberger, A. 
Hamler, M. Trlep, and M. Jesenik, , “Accuracy of 
iron loss calculation in electrical machines by using 
different iron loss models,” J. of Magnetism and 
Magnetic Materials, pp. 269-271, 2003 

[14] Julius Saitz, “Computation of the Core loss in an 
Induction Motor Using Vector Preisach hysteresis 
Model incorporated in Finite Element Analysis,” 
IEEE Transactions on Magnetics, Vol. 36, no. 4, pp. 
769-773, July 2000. 

[15] O. A. Mohamed, S. Liu, and Z. Liu, “Improved FE- 
base PM Synchronous Machine Physical Model with 
Core Losses for Dynamic Simulations,” Proc. Of the 
15th Compumag Conf. on the computation of 
electromagnetic Fields, Vol. 1, pp. 117-118, 
Shenyang, China, June 2005. 

[16] A. Mansouri, H. Trabelsi, and M. Gmiden “Iron 
Losses Calculation of a Surface-Mounted 
Permanent-Magnet Synchronous Motor in Transient 
Finite Element Analysis,” Proc. 4th Inter. Conf. on 
Systems, Signals & Devices, Hammet, Tunisia, 2007. 

[17] A. Mansouri, H. Trabelsi, and A. Masmoudi “Finite 
Element Analysis of a Surface-Mounted Permanent- 
Magnet Synchronous Motor,” Proc. Inter. Conf. 
JTEA2006, Hammet, Tunisia, 2006. 

[18] A. Mansouri, and H. Trabelsi, “Incorporation of 
Vector Preisach Hysteresis Model in Transient Finite 
Element Analysis for a SMPM,” Int. Review of 
Electrical engineering, vol. 2, N. 3, pp. 448-454, 
june 2007. 



Automatic Control and System Engineering Journal, ISSN 1687-4811, Volume 9, Issue I, ICGST LLC, Delaware, USA, June 2009 



hafedh.trabelsi @ enis .rnu.tn 



Renewable Energy and Electric Vehicles Unit, Engineering school of 
Sfax, Tunisia ENIS, Route de Soukra, Km 3.5-BP W, 3038 Sfax Fax : 
+216-74 275 595 

mansouriali2002 @ yahoo ,fr 

^Renewable Energy and Electric Vehicles Unit, Engineering school of 
Sfax, Tunisia ENIS, Route de Soukra, Km 3.5-BP W, 3038 Sfax Fax : 
+216-74 275 595 



Hafedh Trabelsi received the B.S. degree 
from Sfax Engineering School (ENIS), 
University of Sfax, Sfax, Tunisia, in 1989, the 
M.S. degree in the Central School of Lyon, 
France, in 1990, the Ph.D. degree from the 
University of Paris XI Orsay, France, in 1994, 
and the “Habilitation Universitaire” in the 
National Engineering School of Sfax (ENIS), University of Sfax, Sfax, 
Tunisia, in 2008, all in electrical engineering. He is working toward the 
Research Management Ability degree in the field of electrical machine 
design at SES. He joined the Tunisian University, Tunisia, in 1995. He 
moved to the ENIS in 2000. He is currently a Professor of Electrical 
Engineering. He is a member of the Research Unit on Renewable 
Energies and Electric Vehicles of the University of Sfax, where he is 
the chair of the Electric Machine Design Team. He is a member of the 
Organizing Committee of the IEEE International Conference on Signals 
Systems Decision and Information Technology. 




* Ali Mansouri received the B.S. degree in 
electromechanical engineering, the M.S. degree in 
electrical machine analysis and control, the Ph.D. 
degree from National Engineering School in Sfax 
(ENIS), University of Sfax, Sfax, Tunisia, in 2003, 
2009 respectively. He is working toward in the field 
of electrical machine design at ENIS. He joined the 
Tunisian University, Tunisia, in 2003 as an 
Assistant at Gafsa Engineering Institute. 

He is a member of the Research Unit on Renewable Energies and 
Electric Vehicles of the University of Sfax. 



Corresponding Author: Hafedh Trabelsi 
hafedh.trabelsi@yahoo.fr 

Renewable Energy and Electric Vehicles Unit, Engineering school of 
Sfax, Tunisia ENIS, Route de Soukra, Km 3.5-BP W, 3038 Sfax, Fax : 
+216-74 275 595 



41 



Automatic Control and System Engineering Journal, ISSN 1687-4811, Volume 9, Issue I, ICGST LLC, Delaware, USA, June 2009 



42 



Automatic Control and System Engineering Journal, ISSN 1687-4811, Volume 9, Issue I, ICGST LLC, Delaware, USA, June 2009 




i©a®r 



www.icgst.com 






A Qualitative Space Vector PWM Algorithm for a Five-level Neutral Point 

Clamped Inverter 



P. Satish Kumar 1 , J.Amarnath 2 , S.V.L. Narasimham 3 



department of Electrical Engineering, University College of Engineering, 

Osmania University , Hyderabad, A.P, INDIA, 500 047 
2 Department of Electrical Engineering, J.N.T. University, Hyderabad, Andhra Pradesh, INDIA 
3 School of Information Technology, J.N.T. University, Hyderabad, Andhra Pradesh, INDIA 

e-mail: satish_8020@yahoo.co.in 



Abstract 

In this paper a qualitative Space Vector Pulse Width 
Modulation (SVPWM) algorithm for five-level neutral 
point clamped inverter is proposed. With the proposed 
method, the three nearest compound vectors can be 
selected easily and dwelling time of each vector can be 
simply calculated for five-level inverter. The scheme can 
be extended to higher-level inverters. Simulation results 
have been presented for four-level and five-level neutral 
point clamped inverters. The total harmonic distortion 
have been calculated and compared with lower levels. 
The results have been good agreement with the published 
work. 

Keywords: Multi-level inverters, SVPWM, neutral point 
clamped inverter, Induction motor, THD. 

1. Introduction 

The recent industrial applications require high power 
apparatus. The classical Two-level inverters may not be 
useful in high voltage range (above 2 kV) due to limited 
rating of switches. In conventional two-level inverter 
configurations, the harmonic contents reduction of an 
inverter output current is achieved mainly by increasing 
the switching frequency. But the switching frequency is 
restricted by the switching losses in high power and high 
voltage applications [2]. In such applications, multilevel 
inverters have been widely used in recent years for the 
advantage of low harmonic output at low switching 
frequency. At the same time, low blocking voltage in the 
switching devices can be achieved. As a result, multilevel 
power inverter structure has been introduced as an 
alternative in high power and medium voltage situations 
[3]. There are mainly three types of multilevel inverter 
topologies: neutral point clamped inverter, cascaded 
inverter and flying capacitor inverter [4] .The multilevel 
inverters have several advantages like i) Staircase 
waveform quality: Multilevel inverters can generate the 
output voltages with very less distortion. As level of 



inversion increases, the output voltage is nearer to 
sinusoidal, ii) Voltage stress reduction: Multilevel 
inverter can reduce voltage stress on switching devices, 
iii) Input current quality: Multilevel inverters can draw 
input current with less distortion, iv) Reduction in Total 
Harmonic Distortion (THD): Multilevel inverter output 
can have less harmonic content compared to two level 
waveforms operating at the same switching frequency. 
The paper is organized as follows: Section 2 of the paper 
introduces the neutral point clamped multi-level 
inverters. Section 3 illustrates the space vector pulse 
width modulation for multi-level inverters and proposed 
algorithm used to determine the location of the reference 
voltage vector and calculation of switching on-times. 
Section 4 presents simulation results and discussions that 
demonstrate the viability and high performance of the 
proposed algorithm. Finally Section 5 concludes the 
paper. 

2. Neutral Point Clamped Multilevel 
Inverters 

In these inverters, the voltage across semiconductor 
switches are limited by diodes connected to various DC 
levels as such it is called Diode Clamed Multilevel 
inverters. According to the original invention, the 
concept can be extended to any number of levels by 
increasing the number of capacitors addition across 
source dc-bus. Early descriptions of this topology were 
limited to three-levels where two capacitors are 
connected across the dc bus resulting in one additional 
level. The additional level was the neutral point of the dc 
bus, so the terminology Neutral Point Clamped inverter 
was introduced. Because of industrial developments over 
the past several years, the three-level inverter is now used 
extensively in industry applications. Diode clamped 
converter is connected in series-parallel fashion to the 
three phase grid to compensate for distortion in utility 
voltage. The functional diagram of an n-level NPC 
converter is shown in figure 1 . 
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n- 1 




Figure 1 . Functional diagram of n-level NPC Inverter 



To result a defined output voltage, (n-1) consecutive 
switches of each leg must be in the On-state. All possible 
combinations of the above functional diagram can be 
summarized by the equation (1). 



Since its application in 1980s, three-level NPC inverter is 
most practical and widely studied multilevel inverter 
topology. One application of the multilevel diode- 
clamped inverter is an interface between a high-voltage 
dc transmission line and an ac transmission line. Another 
application would be as a variable speed drive for high- 
power medium- voltage (2.4 kV to 13.8 kV) motor. 

The advantages of NPC inverter are: 

(i) All of the phases share a common dc bus, which 
minimizes the capacitance requirements of the inverter. 
For this reason, a back-to-back topology is not only 
possible but also practical for uses such as a high-voltage 
back-to-back inter-connection or an adjustable speed 
drive. (ii)The capacitors can be pre-charged as a group. 
(iii)Efficiency is high for fundamental frequency 
switching. 



^ S .. =1 with i = {a , b , c } ( 1 ) 



The variables Sy are the control functions of the single- 
pole n-throw switches. These variables define the 
position of switches, so that they have the unity value 
when the i output is connected to j point; other wise thy 
are zero. Referring all of the voltages to the lower DC- 
link voltage level (“0” reference) each output voltage 
consists of contributions by a determined number of 
consecutive capacitors: 



Vo = z s # 2> c 

j= 1 V p =1 



with i = {a,b,c) 



( 2 ) 



When balanced distribution of DC-link voltage among 
the capacitors is assumed: 



T , Vdc & , , , 

V i0 = -> with i = {a,b,c) (3) 

n~l m 

The voltage waveform obtained from n-level inverter is 
shown in figure 2. 
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Figure 2. n-level inverter voltage waveform 




In case of five-level inverters four switches from each 
phase-leg will be ON at any point time to produce 
predetermined output at phases. The possible switching 
combinations will be 125 with 64 redundant states. 



Table -I. General characteristics of Multi-level inverters 
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3 -level 


12 


2 


27 


19 


8 


5 


9 


4-level 


18 


3 


64 


37 


27 


7 


13 


5 -level 


24 


4 


125 


61 


64 


9 


17 



a=number of switches with free wheeling diodes 
b=number of consecutive switches of each leg to be in 
ON-state 

c=number of different voltage states of the inverter 
d=number of unique voltage states of the inverter 
e=number of redundant voltage states of the inverter 
g=phase voltage levels 
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3. Space Vector Pulse Width Modulation for 
Multi-level Inverters 

Various Pulse Width Modulation (PWM) algorithms 
have been studied to control the multilevel inverter 
systems and Space Vector Modulation (SVPWM) 
method is a valid one. The most significant advantages of 
SVPWM are fast dynamic response and wide linear range 
of fundamental voltage compared with the conventional 
PWM. But when it is applied to the diode clamped 
inverter and flying capacitor inverter, the SVPWM 
strategy also has to solve the neutral-point voltage 
unbalance problem [5] [6]. There are three main steps to 
obtain the proper switching states during each sampling 
period for the SVPWM method [7]: 

1) Choose the proper basic vectors. 

2) Calculate the dwelling time of each selected vectors. 

3) Select the proper sequence of the pulse 

It is not easy to implement the steps above directly with 
the reference voltage vector amplitude and phase, 
especially in case of high-level inverters. The most direct 
way of calculating the time to decompose all the vectors 
into real part and imaginary part. Another efficient way is 
to calculate the time according to reference voltage 
vector of each phase [7] [8]. Many papers discussed the 
methods to solve the problem and they mainly focused on 
last two processes because the first one is fairly simple 
for three level inverters. But for higher level inverters, all 
the three processes become complex. In this paper a 
qualitative SVPWM algorithm is proposed to implement 
the above process for high-level inverters. 
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Figure 4. Space Vector diagram of Five-level inverter 



In Multilevel inverters the reference voltage vector can 
be reproduced in the average sense by switching amongst 
the inverter states situated at the vertices, which are in 
closest proximity to it. The duty cycles (ON time for each 
state) will be found by equating volt-seconds of reference 
voltage with nearest three states. 



m = d l V 1 + d 2 V 2 + d 3 V 3 (12) 



The space vector (Vs) constituted by the pole voltages of 
inverter V az , V bz and V cz with 120° phase displacement 
can be defined as: 



Where di, d 2 and d 3 are duty cycles of voltage vectors 
(states) (Vi, V 2 and V 3 ) in inverter SV diagram nearer to 
the reference vector and m is the voltage reference vector. 



v, = V az + exp .[j(2ji / 3] + V cz exp .(j{Ax / 3)] (4) 

The space vector, Vs can be resolved into two rectangular 
components namely V a and Vp as indicated below: 

= v a + jV p (5) 
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Equating volt-seconds along a axis gives: 

V s cos axT s =V dc xT 1 + ( V dc cos 60° ) x T 2 (7) 

Equating volt-seconds along p axis gives: 

V s sin a xT s = (V dc sin 60 0 ) x T 2 (8) 

Solving equation (7) and equation (8) gives the values for 
the on- time periods Ti , T 2 and T 0 are 

T = V s x T s x sin( n 7 3 — a) ^ 

1 V dc x s i n ( ^ / 3 ) 



Calculation of duty cycles: 

The vector states at vertices of each region can be 
identified from space vector diagrams shown in figure 4. 
Consider the space vector diagram of sector 1 of 5 -level 
inverter, shown in figure 5. The reference vector m 3 is 
located in region 2 of 3 level inverter. 1114 and m 5 are 
reference vectors located in region 8 of four-level 
inverter and region 12 in five-level inverter, respectively, 
mxl and mx2 (x=3 or 4 or 5) are projections of reference 
vectors on to zero and sixty degrees axes (angle ‘0’ is 
angle made by reference vector from zero axes i.e., 
starting of sector 1; be noted m 3 , irq and m 5 are having 
different angle ’0’ value). The reference vector can be 
synthesized by sequential switching operation of nearest 
three switching states (vertices of the region in which 
reference vector is located). 

The lengths of new vectors can be calculated as 

mi=m x (cos 0 - sin 0 /V3) 

m 2 =2 x m x sin 0 N3) (13) 

The values of mi and m 2 for reference vector in each 
region can be calculated with equation (13). The duty 
cycles of vertices of reference voltage will be m u m 2 and 




45 




Automatic Control and System Engineering Journal, ISSN 1687-4811, Volume 9, Issue I, ICGST LLC, Delaware, USA, June 2009 



[l-(ml+m 2 )]. For example, with reference to m 3 
(reference vector in region 2), the reference vector can be 
synthesized by switching vectors Vi, V 2 and V 3 . It shall 
be important to note that duty cycle for switching state Vi 
shall be length of the vector joining V 3 and Vi, whereas, 
mi is the projection of reference vector m 3 from origin. 
As such, the corrected duty cycle for switching state Vi 
in present case would be (mi -0.25). The length of vector 
joining V 3 and V 2 is m 2 . As such, corrected duty cycles 
for switching states Vi, V 2 and V 3 would be (mi -0.25), 
m 2 and (0.75-m r m 2 ) respectively. 

The values of mi and m 2 are useful in identifying the 
region where reference vector is located, which is the 
major problem in multilevel inverters. The conditions for 
identifying reference vector location in each region and 
the corrected duty cycles for each of the level of inverter 
are shown in Table-II. Once the region is identified, the 
appropriate switching sequence of a region can be 
identified. The ON time period for each state can be 
calculated with duty cycles obtained: 

Tqn for state 1 = Ts x mi 

Tqn for state2 = Ts x m 2 



Tqn for state3 = Ts x [1- (mi+m 2 )] (14) 





Figure 6. Sector 1 of five-level inverter 
Table-II. Location of reference vector in Five-level inverter 



Region 


Condition for 
location of 
reference vector 


Corrected mi, m 2 and 
m 3 for switching 
states 


10 


0.75<mi<l ; 
m 2 <0.25; 
(1+ m 2 )<l 


m,= m r 0.75; 
m 2 = m 2 ; 
m 3 =l- m i - m 2 ; 


11 


0.5< mi<0.75; 
m 2 <0.25 ; 

(mi+ m 2 )>0.75 


mi=0.75-mi; 
m 2 =0.25- m 2 ; 
m 3 = m,+ m 2 -0.75 


12 


0.5<ml<0.75; 
0.25< m 2 <0.5; 
(mi+ m 2 )<l 


mi=m r 0.5; 
m 2 =m 2 -.25; 

m 3 = 1 - m i - m 2 ; 


13 


0.25< mi<0.5; 
0.25< m 2 <0.5; 
(mi+ m 2 )>0.75 


m 1=0.5-111,; 
m 2 =0.5- m 2 ; 
m 3 = m i + m 2 -0.75; 


14 


0.25< mi<0.5 ; 
0.5< m 2 <0.75; 
(mi+ m 2 )<l 


m i = ni] -0.25; 
m 2 = m 2 -0.5; 
m 3 = 1 - m i - m 2 ; 


15 


mi<0.25; 

0.5< m 2 <0.75; 
(mi+ m 2 )>0.75 


nii=0.25-mi ; 
m 2 =0.75- m 2 ; 
m 3 = m,+ m 2 -0.75 


16 


mi<0.25; 
0.75< m 2 <l; 
(mi+ m 2 )<l 


mp mp 
m 2 = m 2 -0.75; 

m 3 = 1 - m i - m 2 ; 



Table-Ill. Switching sequence for Five-level inverter 



Sector 


Region 


ON Sequence 


1 


1 


222 


322 


332 


333 
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322 


422 


432 


433 
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322 


332 


432 


433 
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332 


432 


442 


443 
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311 


411 


421 


422 
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311 


321 


421 


422 




7 


321 


421 


431 


432 
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321 


331 


431 


432 
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331 


431 


441 


442 




10 


300 


400 


410 


411 



Figure 5. Flowchart of the SVPWM for multi-level inverters 
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Sector 


Region 


ON Sequence 




11 


300 


310 


410 


411 




12 


310 


410 


420 


421 




13 


310 


320 


420 


421 




14 


320 


420 


430 


431 




15 


320 


330 


430 


431 




16 


330 


430 


440 


441 




2 


17 


333 


343 


443 


444 




18 


332 


342 


442 


443 




19 


332 


342 


343 


443 




20 


232 


242 


342 


343 




21 


331 


341 


441 


442 




22 


331 


341 


342 


442 




23 


231 


241 


341 


342 




24 


231 


241 


242 


342 




25 


131 


141 


241 


242 




26 


330 


340 


440 


441 




27 


330 


340 


341 


441 




28 


230 


240 


340 


341 




29 


230 


240 


241 


341 




30 


130 


140 


240 


241 




31 


130 


140 


141 


241 




32 


030 


040 


140 


141 




3 


33 


000 


010 


Oil 


111 




34 


010 


020 


021 


121 




35 


010 


Oil 


021 


121 




36 


Oil 


021 


022 


122 




37 


131 


141 


142 


242 




38 


132 


142 


242 


243 




39 


132 


142 


143 


243 




40 


133 


143 


243 


244 




41 


133 


143 


144 


244 




42 


030 


040 


041 


141 




43 


031 


041 


141 


142 




44 


031 


041 


042 


142 




45 


032 


042 


142 


143 




46 


032 


042 


043 


143 




47 


033 


043 


143 


144 




48 


033 


043 


044 


144 




4 


49 


333 


334 


344 


444 




50 


233 


234 


244 


344 




51 


223 


123 


122 


112 




52 


223 


224 


234 


334 




53 


133 


134 


144 


244 




54 


123 


133 


144 


244 




55 


123 


124 


134 


234 




56 


113 


123 


124 


224 




57 


113 


114 


124 


224 




58 


033 


034 


044 


144 




59 


023 


033 


034 


134 




60 


023 


024 


034 


134 




61 


013 


023 


024 


124 




62 


013 


014 


024 


124 




63 


003 


013 


014 


114 




64 


003 


004 


014 


114 



Sector 


Region 


ON Sequence 




5 


65 


000 


001 


101 


111 




66 


001 


002 


102 


112 




67 


001 


101 


102 


112 




68 


101 


102 


202 


212 




69 


103 


113 


114 


214 




70 


102 


103 


113 


213 




71 


102 


103 


203 


213 




72 


202 


203 


213 


313 




73 


202 


203 


303 


313 




74 


003 


004 


104 


114 




75 


003 


103 


104 


114 




76 


103 


104 


204 


214 




77 


103 


203 


204 


214 




78 


203 


204 


304 


314 




79 


203 


303 


304 


314 




80 


303 


304 


404 


414 




6 


81 


333 


433 


434 


444 




82 


323 


423 


424 


434 




83 


323 


423 


433 


434 




84 


322 


422 


423 


433 




85 


302 


303 


313 


413 




86 


302 


312 


313 


413 




87 


312 


412 


413 


423 




88 


312 


412 


422 


423 




89 


311 


411 


412 


312 




90 


303 


403 


404 


414 




91 


303 


403 


413 


414 




92 


302 


402 


403 


413 




93 


302 


402 


412 


413 




94 


301 


401 


402 


412 




95 


301 


401 


411 


412 




96 


300 


400 


401 


411 



4. Simulation Results and Discussion 

The simulation is carried out for the Two-level, Three- 
level, Four-level and Five-level inverters using a 
qualitative space vector PWM algorithm and the 
simulation results are presented. Figure 7 shows the line 
to line voltages of multi-level inverters. It shows the 
improvement of output voltage with the increase of level 
of inverter. The total harmonic distortion is shown in 
figure 8 and change in THD with level of inverter is 
shown in figure 10.. The stator current, speed and torque 
of induction motor are shown in figure 9. 

Simulation parameters are as follows: 

SVPWM parameters: Modulation index m = 0.8; 

DC Supply voltage = 400V; 

No. of switching intervals =192; 
Sampling time Ts =1.0417 e' 4 sec. 
Induction Motor: 2 HP, 1500W, 400V, 50Hz, 

4pole, 1430rpm, Rs= 1.405 Q; 

R r = 1.395 Q; Ls = 5.839 mH; 

L r = 5.839mH; L m = 0.143 H; 

J= 0.0131 Kg-m 2 ; 
fr = 0.002985N-m-s 
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(a) two-level 
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(b) three-level 





Figure 7. Line-Line voltages of Multi-level inverters. 



Selected signal: 50 cycles 




Time (s) 

Fundamental (50Hz) = 176.9 , THD= 28.34% 




(a) three-level 




Time (s) 

Fundamental (50Hz) = 191.3 , THD= 19.97% 




(b) four-level 

Selected signal: 50 cycles 




■o 






5 



Time 

Fundamental (50Hz;> = 19S.3 „ THI> 15.39% 




(c) five-level 



Figure 8. FFT analysis and THD 
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(a) stator currents 





Figure 9. Stator current, Speed and Torque of Induction motor 




Figure 10. Inverter level Vs THD. 



Table -IV. Simulation results with Full load on Induction motor 



Parameter 


Inverter 




2LI 


3LI 


4LI 


5LI 


No. of spikes / 
Half cycle in 
I - steady state 


20 


6 


6 


5 


Peak to Peak I - 
starting 


88 


79 


85 


84 


Peak to Peak I - 
steady state 


19 


17 


16 


15 


Time taken to 
reach I - steady 
state 


0.16 


0.18 


0.16 


0.14 


Peak to Peak T - 
starting 


62 


54 


59 


61 


T - steady state 
(mean) 


10.20 


10.44 


10.3 


10.45 


Speed ripple 


0.3 


0.2 


0.25 


0.5 


%THD 


44.33 


28.34 


19.97 


15.39 



5. Conclusions 

In this paper a qualitative space vector pulse width 
modulation algorithm has been described and applied to 
three-level, four-level and five-level inverters. Compared 
with conventional methods, this method has the 
advantage of ease implementing, especially for the 
inverters with more levels. From simulation results it is 
observed that the generated voltage spectrum is very 
much improved with increase the level of inverter. The 
total harmonic distortion (THD) is highly reduced as the 
level of the inverter is increases. The input current drawn 
by the induction motor is less distorted as level of 
inverter increases. The results have been presented and 
analysed in this paper. 
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